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Abitract 


Spectral Analysis in Geodetic Science and a method for the solution of 
least-squares collocation in the frequency domain as well as a fast Toeplitz 
inversion algorithm in tho s^jace domain are outlined. 

The mathematical background in s|)ectral analysis as applied to geodetic 
u|?olications is summarized. Ingeodetic problems, we usually have discrete 
and non-periodic data. The resolution (cut-off frequency) of the GEOS-3 
altimeter data has been examined by determining the shortest wavelength 
(corresponding to the cut-off frequency) recoverable. The data from some 18 
profiles were used in this study. The total power (variance) in the sea surface 
topography witli respect to the reference ellipsoid as well as with respect to the 
GEM-9 surface was computed. 

A fast inversion algorithm for matrices of simple and block Toeplitz 
matrices and its application to least-squares collocation is explained. This 
algorithm yields a considerable gain in computer time and storage in comparison 
with conventional least-squares collocation. 

Frequency domain least-squares collocation techniques are also introduced 
and applied to estimating gxmvity anoma.'.ies from GEOS-i altimeter data. The. e 
techniques substantially reduce the computer time and requirements in storage 
associated with the conventional least-squares collocution. Numerical examples 
given in this paper demonstrate the efficiency and speed of these techniques. The 
number of numerical operations required to calculate the signals is proportional 
to N logj N (where N Is the number of observations) rather than N® with the 
fast Toeplitz inversion algorithm or with classical collocation. 
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i. Spectral Analyiii and Appllcationi 
1.1 Introduction 


Spectral analyela can be defined ns the process of cnlculatinK and inter- 
preting a S|)cctrum. It actually yields a deeper understanding of data as well 
as the systems which produced those data. The latter property has been very 
usefbl In geophysical applications. 

Spectral analysis is a good tool for the analysis of time or space series. 
This approach is commonly known as "the frequency domain method". While 
using stxHJtral methods one always has to keep in mind that the fre(|uency domain 
a|)|>roach is not the only way to study time series. Time series can also be 
studied using the time domain appi*oach by iMrametiic models such us auton*- 
gressive (AR) or autoregressive moving average (ARMA) models etc. But the 
frequency domain approach luis grt'Ut advantages over the time domain approach 
as we will explain in this chapter. 

The transformation of a rt*cord to the fretiuency domain is U*rmed its 
"spectrum". The frequency and/or wavenumln'r, which is the reciprocal of 
frequency, is in many respects a mort' significant and more useful variable to 
use than the time (or spucel functions. It is to be noted that transforming a 
data vector, say x( t) , into the frequency domain dues not mean the addition 
of anything new but only rearrangement of the given data according to frequency 
(or wavenumber) instead of according to time sequence. The advantage of rep- 
resentation in the frequency domain comes from the fact that most geophysical 
phenomena are expressed theoretically in frequency de|iendcnt forni. The in- 
dependent variable frequency or wavenumber provides a reliable and unique 
check on the data in the time domain and on the systems which produce those 
data, h) the frequency domain, comiiarisons of different records are referred 
to the same value of a certain parameter. Spectral analysis utilizes the whole 
signal, so no information is left out. For details see Bath (1974). 


1.2 Application Fields of Spectral Analysis 
1. 2, 1 Prediction and Interpolation 

As it is very well known, time scries analysis is mainly concerned with the 
study of the time (or space) variations of physical processes. When the state of 
the process is represented by measurements with one or more components at the 
time point, then we can represent the variations of the process over time by a 
vector of real-valued functions. These hinctions can be expressed as a sum of 
an infinite number of sinusoidal terms, i.e. Fourier series. This is nothing 
else but fitting to the data a trigonometric pol}momial in the least-squares sense 
so that the value of the function can be predicted everywhere along the profile or 
surface. 


i.2.2 Amlviii ot Sample DtUt 


W« rtduoe a oompUoatcd hmotlon to a aorlea of simple trigonometric 
fUnotionc, In other «vorda to Its hindamontal wave ftmotlons which provide a 
different and rather revealing Information for the sample sequence. 

A time function x<t) can be characterised statistically when the mean, 
variance, and covariances between the values of x( t) at different times are 
given. But in the analysis of a finite length of record the spectrum, which is 
the Fourier transform of the auto«*nvariance fbnetion as shown later, is often 
preferable to the autocovariauce function mainly for the following reasons t 

a) Estimotes of the spectrum at nei[dd>orlng frequencies ore generally 
uncorrelated. Thus we can interpret the sample spectrum easier than the 
sample autocovariance function. 

b) As it was mentioned before, in many physical problems the spectrum 
is of direct physical interest to us. 


1. 2, 3 Aid to Computations 

VVe can use spectral analysis as an aid to computations, because some 
calculations can be performed more easily In the spectral domain than in the 
time (or space) domain. The frequency-domain representations arc often sim 
pier to handle computationally. This feature will be applied in Chapter 4 for 
the recovery of gravity anomalies from geoid heights which can be called the 
solution by frequency domain least-squares collocation. 


1. 2, 4 Filtering and Control 

By spectral methods we can compute the wave components of a time (or 
space) series and remove some unwanted components from the estimates. Fil- 
tering represents a weighted average in the time domain and a single multipli- 
cation in the frequency domain. So filtering is easier to apply in the frequency 
domain than in the time domain. 


1. 2. 5 Differentiation and Integration 

By means of differentiation the relative distribution of power can be em- 
phasiz.ed in the high frequency (Mayhan, 1978, p. 5-27), for example: 

X(t) =1^ Xn‘c‘'“^0< 

n as 


- 2 - 


x*(t) >7 (lna‘o)X t‘ ‘ 

a i ^ 

f «i« » 

.y n„„,,('xo' ' 

when* Xr, are complex Fourier coefficients, 'V. 2^?/'^ i'. the fitndainentul 

angular frequency, and I Here we can * that x" ( t) would have 

more hl(d> -frequency content than x(t) becaus** of th< ap|H*a ranee of n*. 

If we ait^ Interested in deemphnsizing the high -fnK{uenc> content of a 
signal, then we can use integration, for example: 

■13 

Xj(t) f x(t)dt - T e-'a^»' 

X,., -Jx,,Od. . f (-J^)“x. 

. ; S' - -e 

provided Ums average (mean) value is zero, in this case :^(t) would have less 
high-frequency content than x ( t ). 

In addition to the applications aitovc, the frequency domain methods have 
been successively applied for the estimation of transfer functions for simulation 
and optimization of the data vector, for generating new physical theories, for 
pattern recognition, for studying periodic solutions to physical problems de- 
scribed by diffcrciitis! , lations, for ai^roxlmnting non-periodic functions, as 
an operational device foi 8 )Iving differential equations, etc. 


1.3 Fourier S'^ries and Fauricr Transforms 


1.3.1 The Standard Form of Fourier Scries 


In order to represent a data vector in the frequency domain we shall study 
Fourier scries and Fourier transforms first. According to the Fourier theorem 
(Bath, 1974 , p. 26) a periodic function x(t) having a fundamental period T 
and satisfying the Dirichlct's conditions can be represented by an infinite Fourier 
series « 

“ 4^+/ ( n,, cos noJot + hr, sin ntu>t) (1.1) 

2 -bI 

where Wo”2n/T is the fundamental angular frequency, a,, and b„ are Fourier 
coefficients. The derivation of cqn. (1.1) is given in Appendix 1. The Dirichlet's 
conditions can be summarized as follows: 

- x( t) is at least sectionally continuous with finite jumps 

- x( t) possesses a finite number of maxima and minima 

- The integral J x(t)dt should be convergent. 
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Ai an oximplv a typical periodic function la ahown In Figure 1.1 below. 


*(t) 



We can determine the Fourier cwCficients of eqn. (1. 1) through the 
following AtepB: 

- To find the coefficients Up multiply eqn. <1. 1) by cos n ccU • n 0» 1,2, . . . 
from both aides and integrate over the ix>rlod T with respect to t . 

- To find the coefficients b„ multiply eqn. (1.1) by slnnuvt, n 1,2, ... 
from both sides and integrate over T with res(M.'ct to t . 


In the integration procedure above the following orthogonality relations 
have been used: 

r . I T/2 for m = n 0 

J sin mccot sin nccot dt * ^ 0 otherwise (1.2. a) 

I* ^ T/2 for m n ^ 0 

J cos mcoot cos neCot dt ■ i T for m n~0 (1,2.1)) 

' 0 o^erwisc 

Jsin mu;bt cos ncoot dt « 0 for all m and n (1.2.c) 


Following the steps above, we have for the Fourier coefficients of eqn. (1. 1) 

a,. » yJ n-0,1,2,... (1.3. a) 

2 r 

bp ^ slnnwotdt, n»l,2,... (1.3.b) 

As can be seen from the expressions above: Sp in the cosine transform of x( t) 
and bp is the sine transform of x(t). It is common to call a,, "co- spectrum" 
and bp "quad-spcctrum" (Bath, 1974, p. 40). 


1.3.2 The Complex Form of the Fourier Series 

We can also express the Fourier series in complex form by substituting 
the following identities 
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COB nust ( O' * *»“■ w ’ ) ^2 

Rlnna'. t (»." u> * - c"* ' “» ’>/< 21) 

(n cqn. (1. 1) to o)>tain 

I . ••• 

x<l) ” ***^' ^ ^ ^ ' **' ” **’ *‘‘ 

K(t) / X f 

- SI. * 


( 1 . 1 ) 


(i.r., 


( I . 


whore, 


X. - f (a - ll> ), n 


- 1 . 0 , 1 ,. 


( 1 . 7 , 


are the complex F<jurh.*r ccM*fflclonl». Thin to the complc.x forn (»f tin* FouiMjt 
B orlcB. If wo subBtltuto eqn. (1.3.a»l)) in oqn. wc> Imvc 


X. 




-I I 


cU 


(1. W) 


TIh* eqn, (1. 8) is called the "frequenc-y domain repix senuilan of xi t) '. The 
rosultinK complex Fourier coefficients ar- orthogonal (Pa)!ouli8, 1(10.')) Huch that 


hi 




x.f - 


j 0 tor n * m 
• Of for n ,u 


(1.9) 


whoiv !•; denoUs the e.xpectniioii over T, and L’(Xy)= Eix( t)j, E(X, ) ~ 0 for n f 
We also define a quantity ot,. as follows: 


X r 


nr, “ Y x(t + T) ft*' d 


1 r 

TJ/'( 


)e 


''U-V 


( 1 . 10 ) 


Equation (1.9) above can be proven as follows; 


0) Efx;, 



Eix(t)l 


e 


-1 i.U)n' 


d t 


(1.11) 


since E(x(t)| is a constant this expression is zero for n ^ o 

b) EfX,*x(t)} =• E{x(t)x*(u)le*"^®^lu (1.12) 

where C (t-u) is the autocovariance of the function x, t ) and the superscript ( • , 
denotes the complex conjugate, if x( t) is peritxlic In the mean square sense, 
which will be explained in Section 1. 3. 6, then C ( ^ ) is also periodic, i.c. 

C ( r ) = c ( T + T ). Then it follows 

E{X„**x(t)| == a , e‘'‘*^' (1.1.3) 


Now we can prove that 



EfXn X*] - j: J 

% 

Et' X„* ♦ x(u)3 du 

T 


= J2U 

P e-'(“-")tCo. du = i ° 

n k 

T . 

it 1 n, for 

n “ k 
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Using s similsr approach we can olso show the complex form of the Fourier aeries 
EC|x(t>- f X, e*"*^*!*} - 0 (1.16) 

r.W.ao 

Let us (Eke the square of the homogeneous equation above in order to see the 
cancellations of its terms. 


E{x(t) x(t)}-f E{X,* fX- x*(t))e 


a*o 




'(“-it 


“ C(0)-T ft- e*''' ‘*’ 0 ^ - 2 Or* e"’"*^ + 1 »r. 

- C(0) -5 ty* « 0 (1.16) 

r» mS 00 


The last equality follows from eqn. (1.30) in section 1.3.5. 


Eqn. (1. 8) can also be derived directly from eqn. (1. 6) via the orthogonality 
relationships for the exponential limcticn defined as 

Multiplying both sides of eqn. (1.6) by and integrating on (O.T) we have 

dt = I e-‘“*‘dt 

Interchange the order of summation and integration 

J x(t) dt JC, e"‘“‘^°‘dt = T X, (1.18) 

and 

Xn «=— I* x(t) e“‘'‘^'dt, n = (1.19) 

T T 

Ti:ore is a unique one-to-one correspondence between x(t) and X„ called a 
Fourier transform pair, usually denoted by 

x(t) <&— » x„ 

So for a given x( t) , there is one and only one set of X„ , and vice versa. The 
set of X„ corresponding to x(t) is called the "spectrum" of x(t). 

Since X„ is complex one needs two graphs to display X„. One is the graph 
of the |Xt| versus n, or frequency which is called the amplitude spectrum and 
the other is the graph of the angle of Xn versus n or frequency which is called 
the phase spectrum. The expressions for amplitude |Xn| and phase angle 
are: , 

|X„| - 4 (a,® + b,V (1.20) 

0 „ - -tan"Mb„/a„) Arg X„ (1.21) 
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From these relations we can write 


X. - lX,| e“''^' (1.22) 

As readers already have noticed, sn far (for example from eqn. (1. 1) or 
(1.6) ) the spectrum of a periodic ftinction x(t) is a discrete distribution made 
up of a finite number of frequencies when observational data are used. 


1.8,3 The Cosinewave Form of the Fourier Series 
Still another form for eqn. (1. 1) can be derived usii^ the identity 
ar COB niCot + b,,8in nccot ' cos (nubt -f tO ) 

In eqn. (1.1) to obtain 

00 

x(t) - -K® + 2 A, COS (na'ot + 0.) (1.23) 

where A„ - (a® + bf )^ is the amplitude of the cosine wave and ~ ■ tan"‘(bn/a„) 
is the phase angle of the wave. 


1.3,4 Average Power and RMS Value 

There are some cases where we are only interested in the average power 
contained in any frequency component which comprises x( t) . So here we speak 
about the power of x( t) . Average power for a certain wavenumber or frequency 
is defined as 2|x„ |® and plotted against wavenumber or frequency. This function 
is an even function of n , and contains no phase information. 

The total average power is expressed as 

Pavg = \ x^(t) dt (1.24) 

and the RMS value of x(t) is defined to be 

(x(t)]»,s - Pavg = (”J^x^t) dt)‘ (1.25) 

By Parseval's theorem (the proof is given in Appendix l.A) the eqn. (1.24) can 
be written g 

Pav^g =(^) + i^(a„^ + b,^)/2 (1.26) 

or identically in terms of amplitudes of cosinewaves or by complex Fourier 
coefficients a 

P.f. = (^) 


(1.27) 


- Xo + z£ |Xr.|* - £ lx„|* ( 1 . 28 ) 

RWl na*ti 

In fome reference textbooks eqn. (1.26) Is expressed In terms of average 
degree powers „ 

where P„* - (a* + b*)/2 - A,f/2 - 2|x»r, n»l, 2,..., and Po-(ao/2)®. 


1,3,5 Mean Square Periodicity 

A function x(t) is periodic in the mean bquaro sense if the autocovariancu 
C(') is periodic (Papoulis, 1966), p, .367), l.e, if C ( t) » c ( t + T) for 
every then wo can transform C ( r) to the frequency domain; 

C(T) “ ^ Of., (1.30) 

n 

where a„ are the complex Fourier coefficients of C ( t ) and defined as 

«-> “ T e-“'‘*^^dT (1.31) 

The power spectrum (also called spectral density) of a ftinction x( t) is the 
Fourier transform of its autocovariance (unction given as 

S(w) = j_„C(T)e‘“^^ d r (1.32. a) 

and the inverse Fourier transform can be written 

C(T) = JL r S(W) e’^^dT (1.32.b) 

2 TT J M»(x« 

If we substitute eqn. (1. 30) in (1. 32. a) we obtain 

3) Ct 

where z(cOn) “ e"^^^dr is the Fourier transform of z( t) = exp(inccb'’^). 

00 

If we combine the exponential terms of z(co„) then we get 

z(cOn) = 1 dT (1.34. a) 

This is the Fourier transform of z ( t) = 1 with the angular frequency (uj - ncoc ) . 

Now denote X = a.' - nofe in order to write 

z(X)=J^le'^^dT (1.34.b) 

The equation above is the classical Fourier transform of z ( t) = 1 with angular 
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frequency X . To evaluate the integral above we will show that 

1 <— -> 2n«6((e) 

i.e. 2n'6(u)) and 1 are the Fourier tranwform pairs, where d(U)i is the f 
D irac delta function (unit area impulse) defined as follows: 

(*“ r® 

6 (t )-0 for t^ 0 , I 6 (t) dt = l, and finally, i 6 (t) x(t) x( 0 ). 

•Lop •'•OD 

For the proof, let us compute the Inverse Fourier transform of 2n • 6 (co). 




I e 


:U)t 


d(C 


= JL f” (2n 6(10)1 du) = 1 

2 J«> 00 


(1.34.C) 


Hence, eqn. (1.34.b) can be written as 

1 ► z(X) * 2n»6(X) ® 2Tt»6((o - na*o) 

So we have proved that 

2 (Wn) = f e“‘"^^d'r -= 2 n 6 (co-nc 0 o) 

If we substitute this equation in (1.33) we obtain 

OP 

S(to) = 2tt 2 6(u5-nwb) 


(1.35. a) 


(1. 35. b) 




Eqn. (1.35. b) shows that the spectral density S(co) of x(t) is a sequence of 
equidistant impulses. 


1.4 Finite Fourier Representations of Arbitrary Functions and Windows 

By a Fourier series we can approximately represent an arbitrary function 
x(t) as a harmonic polynomial of finite XN(t) with finite degree N. We have 
to truncate a Fourier aeries at some finite degree, because in practice we cannot 
carry the series expansion to infinity. But how large need N be chosen to have 
a reasonable approximation to x( t) ? In this section we will attempt to answer 
such a question. We will also examine whether 

^llm x^(t) = x(t) 
for each t in the interval ( 0 , T ) . 

Finite representations of arbitrary functions lead us to the concept of a 
**window", which is described as follows; If x(t) and y(t) are periodic 
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ftmcUont with period T such that 

y(t) • J W(t - T) x(^) dT 

then W(t) it called die window through which y(t) views x(t). Some im- 
portant properties, which a window should possess, can be stated as follows: 

1. A window W(t) should be an even function to give equal weight to the 
values of x< • ) about a given point. 

W(t) - W(-t) 

2. W(t) should integrate to 1. 

f W(t) dt = 1 

"T 

BO that if x( t) ■ c a constant, then we would have y ( t) » c , the same 
constant. 

3. W ( t ) should have its maximum value at t = 0 . 

|W(t)| i W(0) for all t 

4. W(t) should concentrate as much as possible around the origin, say t *• 0. 
In order to clarify die concept of concentration, let us consider the area, 

A, defined as .r/s 

A = J W(t) dt, and Ai = 'T/2, 

such that Ax should ^ as close to A as possible. Thus y(t) should 
reflect the behavior of x(t) in the neighborhood of t. 

1. 4. 1 Dirichlet and Fejer Kernels 

Let us define the discrete Fourier representation of Xn ( t) as 

N 

XN(t)=-^ + £ (a„ cos ncOot -f b,. sin ncoot), u)o = 2n/T (1.36) 

and substitute eqns. (1.3.a-b) in the above equation and interchange the order of 
summation and integration to obtain 

XN(t) = ^ J I 2 x(X)(cos ncuoXcoB noJot-fsin ntCoXsin ntOotl|dX+;^Jx(X)dX (1.37) 
which can be shortened to 

Xn{t) r (0. 5 + cos u)o(X- t) + . , . + cos NWb (X-t)) x(X) dX 

tJq 

or as given by many authors 


(1.38) 


wheru IXi(*) ■ 2 (0.5 + cos z + cos 2z * * cos Nz), z co t (1. 40) 

Eqn. (1.40) is identical to 


I>,(z) - sin ( N+ 0.5) 7./sin (0.5 z) 


(1.41) 


The eqn. (1.41) ciui be shown by multiplying both sides of it by (sin z/2) and 
then substituting eqn. (1.40) for D,( z) . This function is even and iieriodic 
with period 2n . The (Unction Dk,( z ) is calh‘d "Diriohlet's kernel" having the 
following roots and maximum value 


IX (z) 


r(2N + l) for z 0 2 ,^ 
1. 0 for z - 


(1.42) 


From eqn. (1. 39) it can be seen that IX ( z ) is the window through which 
Xn ( t) views x(t) . Dirtehlet's kernel has some disadvantages such as not 
concentration satisfactorily about t -• 0 . and slow convergence etc. 


Parzen(1967. pp. 212-213) sum^ests, us an alternative to eqn. (1.36), the 
Fejer's arithmetic mean x,,(t) defined as 


I5.(t) = (Xo(t) +x^(t) + ... + Xr.(t)l (1.43) 

This converges uniformly to x(t), provided x(t) is a continuous periodic 
(unction. 1(N(t) can be expanded into a finite Fourier series as follows 


XN(t) = ^ ( 1 - -|yr^) (a„ cos nuJot + b,, sin nuJot) 
or with integral representation (proof is given in Appendix l.B) 

X,(t) = Y Jq Fn (2n -^~) x(X) dX 

„ 1 r sintN + l) ^z T^ „ 

Fn(z) N + 1 L siniz J ’ z^'2nt/T 

Fn ( z ) is called the Fejer kernel. 


where 


(1.44) 


(1.45) 


(1.46) 


1.4.2 The Modified Truncated Fourier Series 


A modified expression for the truncated Fourier series is defined by 


XN*(t) = x,(t) ( On cos NoJbt + bN sin Nuiit) 


or with Integral representation 


XN*(t) 



x(X) dX 


(1.47) 


(1.48) 
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where 


cC(z) “ D.( z) > cok z - fin Nz/tan (|>) 


(1.49) 


1»4.3, Tukey Means 

Tukey auggestod the following mean to form a harmonic approximation to 
x(t) (Parzen, 1967, p, 215) 

XN(t) = 0,25 ».*(t-T/2N) + 0.6 V(t) + 0.26 x.*(t + T/2N) (1.50) 

We can find the window corrcaponding to Tukey means by 

IX*(-Y^(X-t)-Jl+0.5M^(A-t)) 

+ 0.25 t) + ^1 j x(X) d\ (1.51) 

= jJ^TN||iI(X-t)lx(X)dA (1.52) 

where Ts(z) - 0.25 IX*(z- ^) + 0.5 IX,*(z) +0.25 I\*(z+^), z»2rft/T (1.53) 

and it is called the "Tukey kernel". 

For comparison purposes, we show in Figures 1.2, 1.3 the Dirichlel, 
Fejer and Tukey kernels. It can be seen that the Tukey kernel provides a 
better harmonic approximation, since it is more concentrated around the 
origin than either the Dirlohlet or the Fejer kernel. 

1.5 Fourier Expansion Using Discrete Data 

because of the following reasons we may have to use the values of the 
function x(t) at the N equispaced points; 

- The determination of the Fourier coefficients of the periodic continuous 
function x( t) , with period T , might be difficult because of the difficulty in 
evaluating the integrals defining the coefficients, 

- x(t) might not be known at all points t in the interval 0 t s T but 
may be known only at equally spaced points. 

Suppose the values of the function x(t) are known at N equispaced points 
fi> t; 3 t...i t^.i defined by 

tn = nAt, n = 0, 1, . . . ,(N“1) (1.54) 

where At * T/N. 

Eqns. (1.3. a-b), which are the Fourier coefficients of eqn. (1.1), are 
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represented in dieorete form, m will be proved later, ae foUowa 



a» 

2 

" ITjSo ^0 oo* (2”Jn/N) , n>0,l,...,NH 

(1.65. a) 

where 

b. 

■ 5 x(J *At) sin(2^Jn/N), n ■ 1,2 Nh 

N j mO 

„ r(N-l)/2 if N odd 
1 N/2 if N even 

(l.SS.b) 

Nh 

(1.56) 


In eqna. (1. 5S.a-b) we let n go up to Nh , because in order to know the 
values of the discrete Fourier coefficients (a„,b.) to n <■* , it suffices to know 
them up to N;< . We can Justify this conclusion by using the following identies: 


cos (kN -f n) ub tj ■ oos(nubtj) 
cos (kN - n) Wotj = cos (ncootj) 
sin (kN + n) u'o tj = sin (n oJotj ) 

-sin (kN - n)ob t] ° sin(nabtj) 

where a*o=2rt/T = 2n/N*/^t and tj « J »At and substituting them in eqns (l.SS. a -b) 
above we can sec the periodicity of the Fourier coefficients at basic interval (-Nh, Nh), 
in other words 

u" _ } -Nh s n s Nh and k “ any integer. 

‘'a “ ”kN+B * 


Moreover, we notice that a„ = a_„ and b^»-b.„, where 0 ^ n^ Nh. 


As can be seen from eqn (1. 1) that a„, b^ n ” 0, 1, 2, 3, . . . , Nh are the 
only distinct Fourier coefficients for N equispaced data points, since from N 
equispaced observations we can only determine N other quantities namely a„, 
bg,n=0, 1,..., Nh . 

If we define frequency as the number of cycles in the data interval T(-l) 
then Nh is known as "Nyqulst frequency or folding frequency”. And from dis- 
crete observations we can only determine the Fourier coefficients up to the Nyquist 
frequency. So the eqn. (1. 1) takes the following form for the discrete expansion: 

x(k*At) = (Sn cos 2nkn/N + b„ sin 2nkn/N), for N odd (1.58) 

x(k*At) = ^ (a#cos2n kn/N+b,. 8in2TTkn/N)+^^o8^k, for N even /1. 59) 

A n ssl it 

where Nh is defined as in eqn. (1.56). 


Now consider the following trigonometric identities 


N-l 

V cos 2^km/N cos2nkr/N = 

k *0 
N-l 

i^XgSin 2’^km/N cos2^kr/N = 
sin 2nkm/N sln2rrkr/N = 

k = o 


I N for m = r = 0 or N/2 
0 for m ^ r 

N/2 for m = r 0 or N/2 
0 for all m and r 
N/2 for m = r / 0 or N/2 
0 otherwise 


(1.60) 
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and nuppoM the values of the Amotion x( t) at N discrete points arc given, 
then from eqni, (1, 58) and (1. 60) we can determine the Fourier coefficients 
a,, and b,, which are defined by cqns. fl.55.a-b). i.e. 

2 ^ 

a„ - ^ X x(k*At) cos 2n|(ii/N, n»0,l Nm 

, “,r° <>•»« 

br, “ w £ x(k*At) Vin 2nkn/N, n»l,2 R, 

W k «0 

It is also possible to write the complex and oosinewave forms of the discrete 
finite Fourier series (corresponding U> eqns. (1.6), (1.16), and (1.23)). Suppose 
the continuous Amotion x(t) with square Integrable derivative is defined on the 
interval (0, T) and evoluated at N points 


k»5t*=k»T/N, k »0,1,...,(N-1) 


then the complex form of the Fourier coefficients are given by 
i (a„-ib„) » ^^j,,x(k*Ati 

and the inverse transform x(k*5t) can be expresses as follows 


Nh 




where Nh is defined in eqn. (1.56). 


for N even 
I for N odd 


(1.62) 


(1.63) 


The expression for the cosinowave form of the liinction x(k*At) defined 
at N equispaced points (corresponding to eqn. (1.23) of continuous Ainction 
x(t) ) is 


X(k-4t) . cos(2nkn/N-(®„), k = 0,l (n-1) 

where A„ = (a„* +b„*)^, = -tan** (b>„). 


(1.64) 


The computation of the complex Fourier coefficients defined by eqn. (1.62) 
or the computation of a data vector from the complex Fourier coefficients defined 
by eqn. (1.63) requires roug^y N* complex multiplications and additions. Since 
the time of computation of the transforms above increases by the square of N 
data points, computationally it becomes inefficient or even impossible for very 
large N. For this reason spectral analysis requiring the discrete Fourier 
transform (DFT) has not been so common in the past. 


In 1924 Runge and Konig discovered a method for efficiently computing DFT. 
This method is known as the fast Fourier transform (FFT) algorithm. However, 
Runge and Konig’s FFT algorithm has gone unnoticed and FFT was ''rediscovered 
by Cooley and Tukey in 1965 (Cooley, 1965). FFT takes advantage of the fact that 
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DFT can be carried out Iteratively, which reaulta In a conalderable aavlnga of 
computattcn time and In an Improvement of round-off errora aaaoclated with 
theae computations. FFT requires roughly NlofsN complex multiplications 
and additions of the classical discrete Fourier transform method. A detailed 
explanation of FFT Is given In Appendix l.D. 


The discrete Fourier transform pair for the test Fourier transform (FFT) 
Is given as follows: 


1x1 


1 r** . .. -•.•nur.A 

^ x(k*At) e 

x(k*At) -V 

■wo 

This transform pair is usually denoted 

x(k»At) 


direct Fourier transform 
Inverse Fourier transform 


(1.65) 

( 1 . 66 ) 


We sec that cqn. (1. 62) can be directly transformed into the frequency domain 
using FFT, but not eqn. (1.63), whore the sum is not from zero to (N-1), in the 
present form. By rearranging cqn. (1.63) we have, 


where 


x(k*At) 




1 S n n „ A 

e 



xi"V2 for n-0 

x|*V2 for n - Nh and N ■ even 

Xp) otherwise 


(1.67) 


Eqn. (1.67) Is equivalent to 
x(k*At) - 

where the superscript ( "' ) 
parentheses. 


< 1 . 68 , 

n wO 

shows the complex conjugate of the function inside the 


Now let the sum go to (N-1) such that 



x,k-6i) -‘e z<;' e’“"**^* e‘ (z‘:V 

ns 0 nssO 

(1.69) 

where 

forOsn^Nn 
\ 0 otherwise 

(1.70) 

Finally we can write, 



x(k*At) = Xi(k*At) + Xx'(k‘At) 

(1.71) 

where 
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which can be calculated FFT, and x* (k*dt> la Just the complex conJUKnte of 

X|(k*dt) . 


Further, we can ahow the rclatlonahip between the complex Fourier rep- 
roaentatlona ot diacrete data and oontlnuoua function, i.e. the integral form. 
The aampled function x(k*At) is expreeaed from eqn. (1.6.1) ae (su|)|)oie N 
la even) 

X(k»At) - Z xfV*®'’ " (1.72) 


If we compute the Fourier ooefneienta of the continuoua function x( t) uaing the 
Integral form cqn. (1.16), wo get 


X,, » L f x(t) (1.73) 

Suppoae wo evaluate x(t) at the pointa (k*At), k « 0, 1, ..., (N-1) using X 
computed by cqn. (1. 73). We then obtain 

x(k'4t) ■ f X„ 

niKfoo 

« E Xn e‘-^ " (1.74) 

ti Wm 


By the periodicity of the complex exponential we can write 

^ A 

x(k*At) - Z e“’'’'''^(X,. +(X„-N + Xn.,.) +(Xr.-;,N + X^;;n) + ... I (1.75) 
Equating the coefficients of exp (2nikn/N) of cH|n. (1.72) and (1.75), we gi>t 

xt^ - X„ + (X„. ,N + X..4 ,n) (1. 76) 

From now on will be denoted simply na X„ . 

It is very important to note that the coefficient for the nth frequency of the 
function defined at N points ia the sum of the coefficients of the continuous 
function at the (n, n<fN, n-N, n+2N, n-2N,... ) frequencies. The frequencies 
(n±N, ni;2N, . . . ) are called the "ALIASES" of the nth frequency (pMUcr, 1976, 
p. 119). We already know that the frequency ( N/2) is called the Nyquist 
frequency (see p. 114). The aliases of an observed frequency are those frequem* ies 
V it/ch are obtained by adding or subtracting integer multiples of twice the Nyquist 
frequency. 


1.6 Breakdown of Variance from Sampled Data (at N Points) and Perlodogram 
1. 6. 1 Breakdown of Variance 


The total average power (or the variance) of the function x(k *At) is defined 
as follows 
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(1.77) 


P»»** 


1 

N 


n-i 

E U(k»At)|* 

ll mO 


or In terms of Pourter ooofflolente, i.e. in the frequency domain (from the 
Parseval'B teeorem) 


P«** 


(So/2)* ♦ f X (a,*-*-b/) ♦ 

^li tei • 

Xo*+ jxj* ♦ Pk,^®, where 


(1.78) 



(«^/+»\*)/2 if N odd 
If N even 


If we substitute ( a® b*) » A/ « where A« is the amplitude of the cosinewave 
as defined by oqn. (1.04), we get 


P«®* “ Ao* + 


i X 

nwl 


A® + 


P.. 


(1.79) 


The degree power of this fttnotion is given by 


P 


a 

n 


(a„® + b,®)/2-A,®/2-2|xJ® for 
fto /4 for 

<«v* + bN^^/2 for 

On’ /4 for 


n ® lf2,..«f (N h ” 1) 
n “ 0 

n ■ Nm and N odd^^* 
n = K. and N even 


so finally it is possible to write 



(1.81) 


1.6,2 Perlodogram 


The periodogram is a function of frequency (n) and used very frequently 
in practice. As we will show later it is the discrete Fourier transform (DIT) of 
the autocovariancc function of the data vector In question. Moreover, it is directly 
related to the degree power defined in the previous section. 


The Fourier coefficients of eqns. (1.58-59) can be considered as regression 
coefficients (Fuller, 197C, p. 276). Then by the standard regression analysis we 
can partition the total variance (N*P,v®) of the N observations as follows: 


In(») 


(N*(a® + b®)/2 
IN* Oo®/4 
N* a®/4 
N*(a^+bN®)/2 


for n — 1, 2, . . . , (Nn** 1) 
for n = 0 

for n = Nh and N even 
for n » Nh and N odd 


(1.82) 


where the Fourier coefficients are given by eqns. (1. 55.a-b). Notice that eqn. 
(1.82) is equivalent to lN(n) ~ N* P®. 


( n) is usually known as the ’’paridogram". if we assume that x( t) is a 
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•cc|uenoe tit normiil, Independent random vartablea, thcMi Uu> I'ourivr 

coeffiolonta (a,,b.), being linear combinatiuna of x(t), will be normally 
dlatrlbuted. Thcae oocfflrlenta are lndi«pcndent, becauae the line and coalne 
ftinctlona are orthogonal in the baaio Interval. Aa a reault we can write 

Mn)/c* - Xg® for n ’ 1,2 (N„-l) 

l.(H,)/a® ~ Xi* for N even . ^ 

l.( 0 )/a' 

MNh)/o® -X,® for N odd 


l 4 ?t ua aubatltutc cqna. (l,S6.a-b) In eqn. (1.82) in order to gf>t an alternate 
expreaalon for the periodognim 

MO) •5y^J^|x(k*a)|® (l.Hl) 

M») "■ -jf II ^J]^x(k*At) co8 2nkn/Nj +| ^ x( k*At ) sin 2 nku /Nj j(l. h 5) 
If N even, then 1 ,(Nh) la given by ” 


( Nm ) 


1 

It 


L’ x(k*At) 


K » 0 


C08 



(l.Hfi) 


If we eonaider the following trigonometric identities 

N N 

x(k»At) COB 2nkn/N “ |x(!:*At ) - |i ) cob 2nkn/N (1.H7) 

fi N 

J]x(k»At) Bln 2nkn/N “ (x( k*At) - ) Bin 2nkn/N (1.«K) 

K ■ 1 K a 1 

th(*n eqn. (1. 85) can be written as 

Mn) I*-’*)** 2’tkn/N J +[ J’,^(x(k*Al)-ji)8in2''kn N | 

“ Slx(k‘At)-#i|(x(j*At)-p) COB 2Mn(k-j)/N| 

W '•kts J J -!l '' 

Now let k - J ■» p, to yield 

TT E 5 (x(J*At)- p|(x(j»At + p‘At) -pj COB 2"np/N 
where 1 j+p s N such that 

H, for p » 0 and h for p‘ 0 

substitute these in eqn. (1.90) in order to get 

N-i f2VJ‘l I 

Mn) = 2 J cos 2nnp/N{TT2J|x(J*At)-(i)(x(j*At+ |p|'At)-plf 


(1.8U) 


(1.90) 


Mn) - 2 V y(p) cos 2nnp/N, n«l,2,...,NH 


(1.91) 
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1n(0) - 

where \ ■ [x(At) + x(2At) •»>...■•• x(N/^t)l/N li the mean value of the 
ftmctioD. 


It can be proven (Fuller, 1976, p. 279) that 


where 

and 


Urn E[l 4 (W)l ■ 4rr f(W), CO ^ 0 
Urn E { 1,(0) - 2NM®} - 4nf(0) 


(1.92) 


Having found the Fourier coefflcientB, tctal average power, and degree 
power, now we can compute the uncertainties of these quantities. If the discrete 
observations x(k*At) are independent and normally distributed N(tti<^^) 
Fourier coefficients (an,b„) defined as in eqn. (1.61), i.e. 


2 N-l 

E x(k*At) cos 2Tikn/N 

0 

b„ = "If x(k*At) sin 2nkn/N 

k » 0 


(1.93) 


then, as derived in Appendix l.C 

var(a„) = E {a„ - E(ri„)f 


" N .ICqCos® 2nkn/N 


a N_i 


(1.94) 


using the eqn. (1.60) we find 


var (a„) 


2a^/N for n ^ 0 

4a^/N for n * 0 

4o®/N for n » Nh and N even 


(1.95) 


similarly. 


var (bn) = E fb„ - E (b„) f = { for n ^ 0, or 


n = Nh and N even 


Finally, we can compute the variance of the average degree power, P„ , which is 
a function of the Fourier coefficients as given in eqn. (1.80). 

var(Pn) = ff®/N (1.97) 
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1« 7 Nott-Periodlc FunotloM and the Fourier Traniform 


1.7.1 Mathemattcal DeftaiUop 


Ui previous seotioiis we dealt with periodic functions exclusively. Now we 
wish to find out if it is possft>le to express non-periodic functions as a sum of 
fundamental sinusoidal functions. The answer is positi /e, but the *'sum" is 
uncountably infinite (a continuous sum or integral) and each sinusoid has essen- 
tially zero amplitude. However the sum of all these infinitesimal sinusoids will 
produce a non-zero, non-periodic signal. 

To present the decomposition of a non-tx>riodic x( t), (Mayhan, 1978, pp. 
5-30) we may start with a periodic function Xr ( t) and allow the period to extend 
to infinity so that a non-periodic x( t) results. If we consider the periodic 
function shown in Figure 1.4 



from eqn. (1.6) we have 


00 

:r(t) = Z X, 


1 r-.UJo' 


with complex Fourier coefficients, 

1 r 

X„ = yJ-t/ 2 dt 


(1.98) 


(1.99) 


As T increases the number of important frequencies in the expansion of 
Xt ( t) will also increase, but each having smaller and smaller amplitude. If we 
define frequency as the number of cycles in a specified interval T (period) then 
we have for the ftmdamcntal frequencies fo and Wo 

fo=l/T, U)o=2n/T (1.100) 

then, we have for frequencies 

f = n fo ~ 1, 2, 3, 4, . . . for T » 1 

f = nfo~i, 1, f, 2,... forT = 2 n=-l,2,... 
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( 1 . 101 ) 


Thus, as T « • we have xr(t) -* x(t) and 

llm “ w, -«*><' n «f * 

f -sot) I 

80 a general non-periodic x(t) will result. 


The harmonic content of x( t) will consist of all angular frequencies to, and 
a continuous sum, i.e.in integral, will replace the discrete sum. Since u'o is 
the fundamental angular frequency, in other words the frequency spacing, the 
spacing approaches a small infinitesimal amount denoted by d(o as T 


da* 


lim 

T 


2 n 

T 


( 1 . 102 ) 


This being the case, the amplitude at any frequency will tend to zero. But the 
expansion of xr (t) will not be zero, since we deal here with an uncountably 
infinite sum of quantities. Now using the relations above we can write the 
complex Fourier coefficients corresponding to eqn. (1. 16) as 

doj P"* 

Cn = J _^x(t) e dt (1.103) 


X„ 


It is common to define 


(1. 104) 


to obtain for eqn. (1. 103) 

X„ = 4^ X(w) (1.105) 

or in case of discrete data with N observations (ur sampled values) 

where At is the data (or sampling) interval; and as T -* " , so does N. If 
we substitute eqn. (1. 105) in the Fourier (complex) expansion of Xr ( t) , then we 
end up with a continuous sum for the expansion of x( t) defined as follows 

x(t) = -^ X(W) e‘^‘ dco (1.107) 

Eqns. (1.104) and (1. 107) are the Fourier transform pair (also known as Fourier 
integrals) for the non-periodic function x(t) and usually denoted as 


x(t) X(03) 

There is a one-to-one correspondence between x(t) and X(co). It is also 
common practice to use the standard frequency instead of angular frequency. So 
if we substiwte 2nf for u we have from eqns. (1. 104) and (1. 107), 
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(1.108) 


X(f) - x(t) 

f* 00 • STT f t 

x(t) *= 1 X(f) e' df (1.109) 


In the cate of periodic functions we had related the complex Fourier coefficients 
Xn to the amplitudes of the sinusoidal function. For a non-periodic function x( t), 
it is possible to interpret X( tc) as the amplitude density. That is to sny, the 
amplitude of the sinusoid is zero at any one frequency a', however when it is 
considered over any finite frequency Interval dec (or Acu), the amplitude is not 
zero but equal to 


dctf 

2n 


X(W) 


Thus, the harmonic content of x( t) can be considered as distributed over all 
frequencies with a density X( co). 


For non-periodic functions, the Fourier coefficients corresponding to (1.08) 
and (1.09) are given below 


a(f) = J ^ x(t) cos (2Trft) dt 
b(f) = I x(t) sln(2Tift)dt 


( 1 . 110 ) 

( 1 . 111 ) 


The equations above can be derived from eqns. (1.3.a-b) of a periodic function 
as follows: 

H(f) = Urn™ a„ = r x( t) cos 2nntA’ dt 
substitute f = n/T , to obtain 


a(f) “ r x( t) cos ( 2nft) dt 

J.oo 


( 1 . 112 ) 


Similarly b(f ) can be derived. It can also be seen from eqns. (1. 108) and 
(1. 110) that 

X(f) = a( f) - i b(f) (1.113) 


If the function is known only at discrete points, then the relation between 
the complex Fourier coefficients of a periodic function and non-periodic function 
can be drived through eqns. (1.55.a-b), namely 

N-l 

N 

a(n) = lim — a„ = llm Y> x(k*At) cos 2rrkn/N (1.114. a) 

^-^00 « N — *00 kssO 

&(n) = lim b„ = llm ^ x(k«At) sln2nkn/N (1.114.b) 

N-.00 ^ N-.00 k 


We know that a„ is an amplitude, i/t is the frequency increment, then, 
( T »a„) is the amplitude per frequency increment of amplitude density, and 
finally, a(f) is half the amplitude density at frequency '*f" (Rayner, 1971, p. 
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56). Rayner lays "The relatlcmiblp between amplitude denalty and amplitude la 
analogous to that between probability denalty and probability. " 


In section 1.3.2 we e;q>anded a periodic Ibnetlon x(t) with period T into 
the complex Fourier series and proved that the complex Fourier coefficients Xn 
ape orthogonal, l.e. r a « 

EIX.X.1 ^ for n-m 

where X„ and a„ are defined by eqns. (1.8) and (1.10) xfispectlvely. 

If the fiinction is^not periodic, to be denoted by x( t) , then neither is its 
autocovariance, say C(t), Hence for t / 0 we always have 6(f) < C(0). 

If for some f = fi, we have 6 (t^)=C( 0), then C( r) is certainly periodic. 
Suppose a non-periodic C ( t ) is given. For a fixed t shown in Figure 1. 5 we 
can expand it into the Fourier series in the Interval (-P/2, P/2), where P is 
the given interval length (Papoulis, 1965). 




Figure 1.5. A Non-Periodic Function. 


The coefficients of the expansion would depend on t in the following way, 

|t|cp/2 (1.115) 

Pn{t) = C(T-t) e-‘'''^o‘dT (1.116) 


The non-periodic function x(t) can be expanded into a Fourier series in the 
Interval (-P/2.P/2) as follows 

Mt) 

We can show that x(t) = x(t) in the defined Interval. That Is to say 

E{|x(t) - x(t)|®] = 0 for |t|-^P/2 (l.llS.a) 

where E denotes the expectation. Eqn. (1. 118. a) is equivalent to 

Ef|x(t)|®]- Efx(t) ^t)}- Efx*(t) x(t)] + E{lx(t)|®l (l.llS.b) 
and we know that 
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for |t| < P/2 


E[|x(t)l*| - El|x(t)|»l - 6 ( 0 ) 

We can alao ahow that 

E(x(t) x*(t)l - E(x*(t) x(t)) ■ 6(0) for |tl <;■ P/2 

For the proof: first multiply built sides of eqn. (1. 117) by x*( t), then take the 
expectation « 

E(x(t) x*(t)} = 51 EfX„ x*(t)} e‘'‘Wo‘ (1.119) 

n ■•CD 

and use eqn. (1. 12) to obtain 

E{X,x*(t)} «-~|E{x(f) x*(t))e"**‘»^^dT— C(T-t)e’‘*Wrt'^dT = fl,(t) (1,120) 
finally substitute eqn. (1. 120) in (1. 119) to conclude the proof. 

E{x(t) x*(t)} = 2 d„(t) e'^^^o^dt = C(0) (1.121) 


Similarly we can show that 

E fx*(t) x(t)} = C(0) (1.122) 

so we see that each term of eqn. (1. 118. b) is equal to ^ (0 ). Thus the proof is 
completed. 

As to the non orthogonality of the complex Fourier coefficients of non- 
periodic functions; we will show that they are approximately orthogonal in a 
finite domain. 

First, let us compute the expectation of the products of complex Fourier 
coefficients. [ Xn X«*) • We can define this product as follows: 

x„ X.*= if X„ x*(t) e‘*‘*^o‘dt ( 1 . 123 ) 

now use eqn. (1. 120) to obtain 

E{X„ X,*} = i f 8 „(t) e"‘^'o' dt (1.124) 

We see that eqn. (1.124) is not zero for n 7^ m. i.e. The complex Fourier coeffi- 
cients are no longer orthogonal. If d(T) were periodic with period T, then we 
could easily show that 

0n(t) = «n exp (-incoot) 

and we could write 

VI 
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Since 6 (t) is not periodic the ooeffio lent! 3^ are not orthogonal. However, 
in eqn. (1. 120) if we replace the integration limita by i; • , then we commit 
an error 6( t, P) which la small for large P , so we can write 

P*/lft(t) — f &(T -t) e""*‘*^^dT + 6(t,P) ■ e“‘’‘Wo‘s(nWo) + 6(t,P) (1.125) 

V.. <9 

Substitute eqn. (1.125) in (1.124) to obtain 

P«E{X,X,*} - J e‘<“"^‘^*dt + 6(t,P) e‘*‘^'dt (1.126) 

From the equation above we can write 


Urn P*E{X„X?} - 

^-40P " 


f 8 (n Wo) tor n « m 
I 0 otherwise 


(1.127) 


We see that E[ |Xn|®l tends to zeroes (1/P) for P- ®, whereas the 
autocovariance E[Xr X.*} tends to zero faster, i.e. the correlation coefficients 
of X„ and X, tend to zero also. So we can say that the coefficients Xn's of 
eqn. (1.117), for large P, are approximately uncorrelated. 


1.7.2 Properties of the Fourier Transform 

The Fourier transform has some very important properties such as line- 
arity, scaling, time shifting, modulation, symmetry, differentiation, convolution, 
etc. , which are well explained in every text book in this topic (e.g. see Bath, 
1974, 11. 42-48). As a result of these properties, some important relations are 
given in Table 1. 1 . 


Table 1. 1. Summary of Fourier Transforms. 


Function 

x(t) 

X(t) 

aiXi(t) +aaXa(t) 

x(at) 

X(t-to) 

exp (1 Wot) x(t) 

X(t)*C08Wot 

Xi(t)*Xg(t) 

Xi(t)*Xa(t) 

p"x(t) [n“ derivative x( t)J 
x(-t) 


Fourier Transform 

X(w) 

2n x(-w) 

aiXx(w) + aaXa(to) 

(l/|al) X(w/a) 
exp (-Iwto) X(U5) 

X(w- Wb) 

( X(w- Wo) + X( w+ Wo))/2 

Xx(w)*3^(w)/2n [convolution] 
Xt(w)*Xa(w) 

(iw)"X(w) 

X(-w) 
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In the following ■eotlone Qie convolution property will bo applied. The 
convolution of aibitrary ftinctione x(t) and y(t) ia defined aa follows: 

x(t) *y(t) - x(T)y(t-T)dT 

- f x(t-T) y(T) dT 

Wo can apply the convolution both in the time domain and in the frequency domain. 

1. 7. 2. 1 Convolution in the Time Domain 

If x(t) X«c) 

y(t) Y(u’) 

then x(t)*y(t) X(aj) Y(U3) 

The above transform pair can be proven by the following procedure: 

A 

z(t) = x(t)*y(t) «= I x(T)y(t-T)dT 
The Fourier transform of z( t) is 

Z(w) = F(x(t)*y(t)J = (“ e-‘<0‘f* x(T) y(t-T) dTdt 

00 v— 00 

H tx> aQD 

= I x(T)dT j y(t-T)e“*^‘ dt 
J_,n 

Now let X * t - T , dX = d t and with t held constant 

Z(M) = J.„XC') dX 

= H *(T) e-‘“’dT J_^y(X)e-'«^ dX 
Finally we can write 

Z(uj) = F(x(t)*y(t)l - X(w)Y(u;) 

1. 7. 2.2 Convolution in the Frequency Domain 

If x(t) X(OJ) 

y(t) Y«c) 

then x(t)y(t) l/2n X(to)*Y(w) 
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■ ■■ 


We CM prove this trenaform pair by defining x(t) > x( t) y(t) and taking the 
Fourter tranaform of s ( t) 


Z(u)) “ Flx(t)y(t)I - f x(t) y(t) e'“*’‘ dt 
For x(t) aubatltute the following Inverae Fourier tranaform 


to obtain 


1 f* iXt 

*<*> ’T;rJ..x<^)« 


z(“) - ““ 


Since the second integral Is the Fourier transform of y ( t) , (u)-X) being 
the angular frequency, we have 

Z(te) = Flx(t)y(t)l J_“„X(X) Y(uj-X) dX 

or equivalently I 

Z(o)) - X(W) * Y(W) 


1, 7, 3 The Finite Record Length and Windows 
1» 7. 3. 1 Time Domain Windows 


In practice we cannot record the results of an experiment of infinite length, 

so Uie use of limited record lengths are inevitable in spectral analysis of obser> 

vational scries. This forces us to truneutc the record at some length, say T , 

such that , ^ - .1 

_ r x(t) for |t| T/2 

"to otherwise 


y(t) 


(1.128) 


where x(t) is a non-periodic function defined for (-®,®) but recorded only 
in the interval (-T/2, T/2), this record being denoted by y(t). Actually this 
truncation is nothing else but applying the rectangular (box-car) window. Then 

we can write ltl 'T/2 

y(t) = X(t)u,(t). where ,l.J29) 


thus, the Fourier transform of y(t) is 


Y(X) = J „y( t) e"‘* dt = J_^^^y(t) e- dt 

Y(X) = J”^x(t) Ur(t) e-‘^‘ dt 


or equivalently 


(1.130. a) 

(l.iao.b) 


If we substitute eqn. (1.107) for x(t) above and evaluate the integral, in 
other words, use the convolution technique, then we have 

Y(X) = -~X(X)*U(X) =-^J“x(w)UT(to-X)do.' (1.131) 


- 28 - 


wherei UT(a»-X) is the Fourier transform of ur(t) defined by 


Thereforoi 


Ut(u)) 

Y(X) 


(Sin u) T/2 ) /( (u/2 ) 


1 

T7T 


1 r* tin IM-A»T/a)l 
(li-A /i- 


(W-X)/2 


(1.132) 

(1.133) 


In the limit as T » the rl^t-hand side of oqn. (1. 132) reduces to X(a*) . 
Our main goal is to determine the true spectrum X((u) » but because of limited 
data we can only estimate Y (m) . So we must design a window, say w ( t ) , in 
such a way that Y (te) will bo very close to X((U) . Even with the application 
of the most optimum windows wo will get a distorted spectrum. The main task, 
then, is to keep the spectral distortion close to a minimum. 


Windows can be applied either in the time domain or in the frequency 
domain. In the time domain; 


y(t) » y(t)»w(t) = x(t) (ur(t)*w(t) 1 = x(t)‘w(t) 

In the frequency domain ; 

Y(w) = X(u))* W(o)) « V(o;)*W(a}) 

That Is to say, eitlier we compute Y(u)) and convolve It with W(u;), or instead 
we multiply the data vector by a window vector in the time domain and compute 
the Fourier transform. 


We can state some of the properties, which a spectral window corresponding 
to the applied time window must possess: 

- small or insignificant side-lobes, i.e. smooth time window without 
sharp corners 

- high concentration at the main (central) lobe 

- symmetry with respect to the y-axis, in other words, the window 
function must be even. 


Designing an optimum window for a specific purpose has been a major 
challenge to scientists in this subject. 

Bath (1974, p. 157) classifies windows as follows; 

(1) Trigonometric windows, which use trigonometric functions of time 

(2) Power windows 

(3) Exponential windows 

Some of the most commonly used windows will be explained below. 
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(1) Rect»nguUir<Box»car or Barlettt Window 

I for |t|*T/2 


{ 


0 othonvtM 


(1.134. a) 


(1.134.b) 


w(t) ' 

with speotnil reprefontetion 

W(uj) « (Sin coT/2)/( w/2) ■ T alnc (wT/2n) 

where slno (x) >■ ■in nx/( nx) ia how the elnc function is defined. A rectangular 
window for a Discrete Fourier Transform (DFT) is defined as (Harris, 1978, p. 

58) 

^ ^ f 1 for n»0,l (N-1) ,, 

«'"> ’lo Otherwuc (l.m.a) 

then the spectral window for the DFT window is given by 

“ ' ”Bln^2)^ e-*^-^, 8=0,1 (N-1) (1.135.b) 


Since die sinc-fUnction has large side-lobes, this window is not very suitable 
for a true representation of a spectrum. In addition to this, the magnitude of these 
side-lobes decreases slowly and half of them are negative causing displeasing 
results, because the power is positive by definition. 

Notice that a large T will lead to more details, i.c, to better resolution 
in the computed spcctnim than a small T . But a small T will lead to better 
stability and reliability of the computed spectral estimates, since for small T's 
the spectral smoothing extends over a larger frequency Interval. 

In order to eliminate the effect of the rectangular window and the effect 
of noise, Bath (1974, p. 180) suggests to use smoothing in the frequency domain. 



where X„ are complex raw Fourier coefficients, whereas x' are smoothed 
complex Fourier coefficients. The above expressions are equivalent to applying 
the Hanning and Hamming lag windows in time domain respectively. 


(2^ Triangular iFeler) Window. Afft> 

1 -2|t|/T for |t| T/2 

0 otherwisv 


4r 


A-(t) ■ {■'T'".'":.!'.:..."* 

with tho ipcctml representation 

. , T sin*(u} T , 8 

**'‘*'> 2 ( CO 2 

The triangular window fora dfT la defined as (Harris, 1978, p. S9i 

•uni ^ for n» 0,l (N^2) 

U(N-n) for n - (N'2)....,(N-1) 

with the correaponding spectral window 

" N • * "•' 

This spectral window has no ne'^ativc side-lubes, but it has large main lobes, 
which do not have high concentration. 


(l.ilifl.b) 


(1.137.u> 


(3) Hannina (Tukey or Cosine) Window 
. tl +COS 2nt/T) /2 


for |t|‘^ T/2 
otherwise 


(1. 138.a) 


v.ti. I . ctral representation 

aincc(T/2) T r sln(U)(T/2) -«-n ) sln(oc(T/2)- 
,tt!/2) +"41 u!(T/2)+n + ct,>(T/2y-'' J 

•“8inc-^ + ^[sinc(co(T/2n)+l) +8inc(u,'(T/2^-!)-l)'j (l.l38.b) 


or as a function of W (co) of rectangular window 


Hj(cc) = 1 W(a*-“-)+ | W(u?)+ i W(to + ^) 

The Hanning window for a DFT is defined as (Harris, 1978, p. 60) 

fi^(l -cos2nn/N) for n ” 0,1, . . . ,(N-1) 


h,(n) = I 


otherwise 


with the corresponding spectral window 


(1. 138. c) 


(1. 139. a) 


Hi(0) =iw(e) + ^(W(0-^) + w(e + ^)l (i.i39.b) 

where W(0) is defined in eqn. (1.135.b). Thus, we sec that Hi( 0 !) is a sum 
of three sine -functions displaced relative to each other and the side-lobes from 
the three sinc-functions cancel each other to a large extent, which is pleasing. 


- 31 - 


<4^ Himmlag Window 

hi(t) - M0.64+0.46coi2nt/T) for |t| -: T/2 
I 0 otherwise 

with the epeotnil representation 


(1.140.S) 


-0.54Tslno|^ ♦0.23Tjslne( 1) +slnc( - 1) j (1.140.b^ 
or In terms of W(u;) of the rectangular window 

Ila(u>) '0.54 W(«) + 0.23( W(u; +(2r/T)) + \V(a)-(2’i/T)) ) (1. 140. c> 


The Hamming wludo^v fora DFT Is defined as (Harris, 1978, p. 62) 

f0,54-0.46cos2nn/N for n ■ 0, 1, . . . ,(N-1) 

I 0 otherwise 

with the corresponding spectral window 


(1.141. a) 


Ha(®)' 0.54 W(e)+0.23(W(fl-(2n/N))+W(e-f(2n/N))], (1.141.b) 

d == 0, 1 , . , , , ( N~l) 

where W(6) is defined in cqn. (1. 135.b). Notice that the Hanning and Hamming 
spectral windows correspond to a weighted average over three consecutive values 
of the spectral fbnctlon of the box-car window. The Hamming window also has 
small side-lobes. 


(5) Cosine-Tapered Rectangular Window 

This window Is flat over most of the dab), but tapers off near the two ends 
of the data. For example, it may consist of cosine window at the ends ( say 10% 
at each end) and a rectangular window in between defined us follows; 


/(1+cos10tU/T)/T for -(T/2) t • -0.4T 
c(t) = J 2/T for ~0.4T<t ' 0.4T 

\ ( 1 + cos lOo t/T)/T for 0 . 4T •' t - T /2 

with the spectral representation 

1 aln cc(T/2> -t-siti 4co(T/^Oi 
C-(W) i: a?(T/2)Il-(W(T/l0n))‘') 

The above window for a DFT is defined as (Rayner, 1971, p. 83) 

rid -cos (lOn n)/N 1 , 0 ' n- N/10 

c(n) =il , N/10 • ir' 9N/10 

U[l“Cos(10rr(N-n))/Nl,9N/l0- n' (N-1) 


(1. 142. a) 


(1.142.b) 


(1.143) 
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An uxU*nilvv ntudy of vorloun window* dcilgm^d for different |iur|)one* I* 
given In Hnrrln (1978). 


1 ♦ 7, 3, 2 Frequency Domnin Wlnttowi 

It !• alto possible to apply windows in the frequency domain and compute 
tht* effects In the time domain. As an example, su|ipost> that X(u;) of a function 
x(t) it tnincated outside (>P, P) yielding X'io;) defined as 


X(w) « X(uJ) W,(u», where W,(u» « 


(l.Hl) 


Since eqn. (1.144) corresponds to a convolution of x(t) and w (t) in the time 
domain, we have 

Xo(t) » x(t) • w„(t) - x( ~ d T (1.14'i) 

Instead of the rectangular truncatior above, if we eliminate the components 
of |a*| > P and favor linearly the low-freciuency com|K>nents, then we hJive 


X^(a') X(u,') • Ap(UJ) 

where (u*) is the triangular pulse deftned by 


(1. 146. a) 


»p(o.) Q otherwise 


(1.146.1)) 


The corresponding equation in the time domain can be exprc‘ssed via the con- 
volution theorem as 


X,(t) == 

^ J.® TT(T/2)(t-r) 


(1.H7) 


1.7.4 power Spectrum Variance Breakdown of Non-Periodic ^^^nctlonH 

The energy content of a non-periodic function x( t) is defined by 

e'* = J*_” x"*(t) dt (1.148) 

or equivalently In the frequency domain by, 

if - ■— J_„X(w) X(-u))da> = — J ^|X(UJ)(‘* dw (1.149) 

If x(t) docs not have a finite energy, then we have to consider the total average 
power P® expressed as follows 

T 

P® - |lm^ ~ x-’( t) d t (1. 150) 

Papoulis (1962, p. 240) classifies functions into the following three classes 
according to the magnitude of the power: 
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(1) P* ■ 0 , obviously oontslns all ftinotlons with finite enem content . 

(2) 0<P^<«, contains all periodic ftinotlons, i.e, P* *,]jrjx„r 
where Xn are complex Fourier coefficients of the function x(t) , 
and of some non-perlodlc ftmctions as well. 

(3) P*-», i.e. x(t) has infinite power. 

Power Spectrum from Sequenced Data 
(aj Using Autocovariance FVinctlon 

For the statistical analysis of dependent, sequenced data, say x(t) , die 
mean variance and the autocovarianoe ftmotion characterize the distribution com- 
pletely. But generally the autocovariance ftinction is difficult to inteipret and 
adjacent estimates are not independent, so the confidence intervals are difficult 
to calculate. On the other hand, the statistical estimates in the frequency domain 
are relatively stable and independent. 


The power spectrum S (to) of a ftmotion x (0) is defined as the Fourier 
transform of its autocovariance function which is expressed as 

= 11m 4: f (1.151) 

' ' r-*® TJt 

or in case of discrete observations 


1 ' 

Then, the power spectrum can be written 

S (it) = J 


(/) = r*A0 (1. 152) 

r = 0, 1 (N-1) 


(1.153) 


Actually S(to) and are the Fourier transform pairs usually shown 
as 

y(^) « S(w) 

The inverse Fourier transform of eqn. (1. 122) givi‘s the autocovariance function 
y((i)) = -^ J_”s(U)) e‘"'l' dili (1.154) 


The power spectrum S (ta) can be directly expressed as a function of 
x(0), (Rayner, 1971, p. 77), 

S(OJ) = llm |X((C)|® (1.155) 

I *400 X 

The relation between eqns. (1. 152) and (1. 155) is known as the Wiener-Khint- 
chine relation. 
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The eqn. (1. 1S2) cannot be applied In ita present form, because y( $ ) 
is not defined outside i> In the case of observed data vector x( lii) . In the 
derivation of eqn. (1. 162) we Implicitly assumed that x( ij>) was zero for 
1 I > T/2, but It Is not necessarily the case. Suj^se 0 Is In the range of 
(•m,m) and let C 0^) be an estimate of the ensemble value of the true auto- 
oovarlanoe ftmotlon V(0), l.e. 

Average (C(^)) - (1.156) 

then, we can make eqn. (1.156) apply to all values of $ through windowing, 
C’(^) - C{$) W(d>) 

where W((|i) is the window fiinctlon. Therefore, eqn. (1. 166) becomes 
Ave(C'(d»)l “ Ave[C((^i)l*W((ii) « y(0)*W((/)) 
and the frequency domain representation 


Ave(S( 0 ))l » S(w) * W(«) (1.157) 

where S(o)) is the spectrum (power) of the true autocovariance function of x( 
(^) is the win^w through which y{$) is viewed, and W(oJ) is the corres- 
ponding spectral window. 

The major steps to compute the power spectrum from a series of N 
equally spaced observaticms can be stated as follows: 


( 1 ) 

( 2 ) 


Remove mean and non -desired effects of low frequencies (trends). 
Compute the unbiased sample autocovariancc function C ( 0) through 


C ((&) 


1 

= x(k + (P), \vhere 


f (ii = r *A((i 
to ^ r s m 


(1. 158. a) 


It is also common tr use the biased sample autocovariance function of the 
following form, 

N-i-|r| 

C(0)=-^ E X(k) X(k+0) (1.158.b) 

N k «B0 


(3) Apply the selected window to get the windowed autocovariance function 


= C(^) W(j6) 

(4) Transform C '((/)) into the frequency domain as follows 

■ 

S'(f„) =^IT^C'(r) cos (2rtf„ r *A0) A0 (1.159) 

where A0 - constant ’"As the lag increment, and As = observational 
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interval, ■ n/<2m* ) wltii n"0, 1,.*., i , (2m*Ai6) ■ the length of 
the baalo Interval. 

Since C^r*h0) and the ooalne ftinotion are even, it ia posaible to write 

I 

S'(fk) ~ 4 r*^il») coa (2nfk r*^|li)A(^ (1.160. a) 

We have to make adjuatmenta for the end data (recall the trapezoidal integration) 
(fbr detalla aee Bath, 1074, pp. 171-173), therefore, 

S'(fk) » 2A4»f(f(0)+2L,C\r*Ail»)ooa-^^-E5li+C|[m*A^)co§ nM (1.160. b) 

r ml m 

with the inverse transform (Otnea, 1072, p. 107) 

C'(r^(|i) » — -i-— rs’(0)+2y’ S*(f„)oo8^^+S*(m)coB Tir) (1.160.C) 

4m nVi m 

The eqn. (1. 160. a) can also be expressed in terms of average variance per 
prequency band (Rayner, 1071, p. 81), by dividing it by the length of the basic 
interval, i.e. (2m*d0) 

S'(fn) “ “ ^ Cf(r»d4>) cos 2nf„r A(ii (1.160.d) 

m rirfo 

The relation between eqns. (1.160. a), the variance density estimate, and (1. 160. d), 
the average per frequency band, are very important and they must be used in the 
right places. 


Still further we have to make adjustments for the spectrum ends thus finally 
we obtain for the power spectrum 

S*(0) = •:— C'(0) 2 C'(r»A((J) C’(m*A((i) 

2m m 2m 

2 rrnr 1 

S’(fn) = ~ C\0)+ — COS— + -c\m*A 0<'n-m 


^'<*”>•■^^'(0)+— (-If C'(r*^»/') + ^^ C'(m*AiP) (1-161) 


(5) Convolve the raw powers above with W( r) , i.e. with the spectral window, 
to get ‘S( r) . As an example, suppose we decide to use Hanning window. Then 
there are two options: (1) we can multiply the data by window function and trans- 
form to the frequency domain, or (2) we can transform the data vector to the 
frequency domain, in other words, a rectangular window is used in the time 
domain, and then airily convolution in the frequency domain, which takes the 
following form for this particular example. 

S(0) = i S’(0) S*(l ) 

5’( r) = i S'(r-l) S’( r) + i S'(r+1) 

§^(m) = i S'(m-l) S'(m) ' ’ ' 
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The above eetimatee are a/erages for bands l/(2m*A^) wide, except 
at r s 0, m where fhey are l/(4m*^^ ) wide centered at cc n r/(m*A0) . 


fbi Using Direct Data Transform 


In order to compute the power spectrum we can also transform the data 
directly. If x(t) is a non-periodic function and only defined in the interval 
(-T/2, T/2) , then we have to use a window, say w(t), to obtain 

dt| 

= ^1 X(f) * W(f)|* (1.163) 

where f - co/2tt , and N => the total number of observations. Now recall the 
following relationship 

X(f) - X(f) * W(f) = a(f) - l6(f) (1.164) 

in order to write 

^(f) = ^[^•(f) + 6^(f)l, |f| < (1.165) 

Since eqn. (1.165) is an even function, for the one-sided spectrum A t* have 

S(f) = ^{a*(f)+fi®(f)l, 0 f (1.166) 

If we convert the power density estimates Into average variance fre- 
quency band, multiplying by l/N, and substitute eqns. (1. 114. a-b) above v;e 
obtain 

^(k) = + (1.167) 

This equation is equivalent to the expression for the degree variance of tJ\.? 
periodic function. But here we have used windowed data and the estimate above 
is an average for a band centered at "k" rather than a discrete estimate at k . 

We see that S(k) contains (N/2) full bands and covers the same range 
of frequencies from zero to (Nyquist frequency) as ^ ( r), which ontains 
m full bands. So we can make S(k) equivalent to lS(r) by summing over 
(N/2m) bands. 

Now we can give the steps for the estimation of the power spectrum from 
a set of observations, x( t) , as follows: 

(1) Remove mean and trends. 
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(2) Apply the window function w(t) , e.g. uie a coalne-tapered rectangular 
window given by eqn. (1.142. a), thus 

ft(t) " x(t)*w(t) (1.168) 

(3) Compute the Fourier coefficients a(f), i^(f) through eqns. (1.110) and 
(1. Ill) or the coefficients A(k), l&(k) by eqns. (1.55.a-b), which are related 
to a(f) and ^(f) through eqns. (1.114.a-b), using F FT in the calculations. 

(4) Calculate the power spectrum d(k) by eqn. (1. 167) and sum over blocks, 
say 5 to 50 (Rayner, 1971, p. 83), to obtain ?(r). If we assume N is even 
such that N/2m 2z-M, then 

^(0) = (a®(k) +b®(k)]/2 

r(8^l)4>t ^ ^ 

0^r-"m (1.169) 

?(m) = (a^(k) +6®(k)l/2 + (-^) 

frequency bands are centered at 

0, l/2m*t, 2/2m*t,..., r/2m*t,..., l/2t 


( c ) Using the Filtering Method 

This method will not be discussed here. In this method data are passed 
through a band-pass filter, squared, and then summed with a final normalization. 
For this method, readers are referred to Otnes (1972, pp. 311-315). 


1. 7.5 The Confidence Intervals 
The chi-square distribution is given by the formula below, 

2 “^° T( u/2 ) 

where r(t^/2) is the gamma function, and u » degrees of freedom. 

It is well known that 

/ DF ♦ sample var lance \ 

\ true variance / 

has X* distribution, DF being degrees of freedom, i. e. in other words, 
u*o*/a^ = x®(U) 

By rearranging ihls equation we can write 


(1. 170) 


(1.171) 
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( u, smaller probability level) ' x"* 1 w, larger probability Itwol) ' 

Now let us see the application of the distribution in case of spectral 
analysis, by going back to the spectrum of the sample variance 

r (f i 

^®(r)-ff(r)« Z , la®(k) +6*(k)l/2 (1.173) 

k»r(^ zn)- t 

Since a(k) and l^(k) are defined from an orthogonal transformation and are 
generally assumed to be normally distributed, 9^( r) of eqn. (1. 173) has 
2(2z 4- 1) degrees of freedom in the interval 0 < r < m , and 2( z -t- 1) degrees 
of freedom at r - 0, m . 

But using a window reduces the degrees of freedom (Tukey, 1967) and 
must be considered in computations, e. g. if we use a cosine-tapered rectan- 
gular window, then 

DF = u - 2 (2z + l)(N-G)/N (1.174) 

where G is the percentage of tapering, and N is the number of total observations. 


1, 7.6 Cross Spectral Analysts and Coherence 
1. 7. 6, 1 C roes Spectral Analysts 
Consider the following linear regression, 

y = ax + b 

where y, a, b, x are functions of "t" . I et y also be a function of all the 
X *s in the region, then we can write 

y(t) » b(t)+ f a(q)*x(t-q) dq (1.175) 

CD 

If x(t) and b(t) are uncorrelated, then by eqn. (1. 175) we obtain 

Cxy(iii) = f a(q) • Cxx(»/»-q) dq (1.176) 

cx> 

with the frequency domain representation 

Sxy(f) = A(f)* Sxx(f) (1.177) 

and 

A(f) = Sxy(f )/Sxx(f) (1.178) 

A ( f ) is known as the response function of the system. 

Recall the definition of the correlation in the time domain, i.e. 
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pv - (2xy/N)/I(J:x*/N)*(EyVN)^l 


(1. 179) 


The correeponding oorrehttlon In the frequency domain la defined as 

- 8*y(f)/((8„(f))*(Syy(f))^l (1.180) 

which is called "the coherence or coherency" and varies between (-1, 1). m the 
equation above, Sv(f) is the Fourier transform of the cross-correlation fbnction 
Cxy(0) and dehned as follows 

8xy(f) » f Cxy((|))* e"**"'* d((), .» «;' f <r 00 (1.181) 

or equivalently, 

Sxy(f) » T (1.182.a) 

and for one-sided spectrum we have 

Sxy(f) =2f Cxy(«i))*e’‘*^^'*'d<i> - Urn I: Y(f)*X*(f) (1.182.b) 

tl.OO f-*00 1 

where Y(f) and X(f) are defined as in eqn. (1.108). 

For discrete and finite data the eqn. (1. 182. a>b) can be written (Rayner, 
1971, p. 82) 

Sxy(f) = .I^(f)»^*(f) = 4(Y’(f) * H(f)l(X’(f) ♦ H(f)l* (1.183) 

-00 <• f <• 00 * 

and for the one-sided cross spectrum we write 

Sxy(f) = .lY(f)‘^*(f) = -|(Y'(f) * H(f)UX’(f) * H(f)l (1.184) 

If cross-covariance density estimates are converted into average cross spectral 
density (CSD) per frequency band, then we have 

iy(k) = -^{•|-[^ay(k)- l|fey(k)][-|5,(k) +l^fex(k)]} 

“ l(ay(k) - i$y(k)Uax(k) + l6x(k)l (1.185) 


The real parts of eqns. (1. 183) and (1. 185) are called cospectra, and 
imaginary parts of those are called quadrature spectra, that is to say 

Sxy(f) = ly(f)“l&y(f) 

or equivalently, (1, 186) 

Sxy ( k ) = Rxy ( k ) — I Qxy ( k ) 
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where Aiy(f) end K^{k) ere ooepeotra, and ^iy(f) end ere qued> 

reture epeotre. 

We cen e>q)reeB the oroee-ooveriencee between x end y (eeeuming means 
ere zero) 

(e ) in continuous oese: 

C*y(#) - lim 4rJ„y(‘)*x(t+v) dt (1.187) 

( b ) in discrete esse: 

1 ^-4' i 

Ciy((/>) - nITTT 4 I ,y(‘)*X(t+<y), (1.188) 

Notice that eqns. (1. 187) and (1. 188) are not event therefore» we break them into 
even and odd series defined as follows: 

(1.189) 

(I.IOQ) 

Eqn. (1. 189) is transformed into frequency domain through eqn. (1. 160) and (1. 162) 
to get S'ty^(r) which is known as "the cospectrum of y(t) and x(t)", and de- 
fined by , 

S’,yj(r) = 2 A(() (Clyj(0)+ 2^^^cly( r) coB-^~ + Clyj(m) cos rr r] (1.191. a) 

The power of the eqn. (1.190) is eomputed similarly, 

■ -i 

2 

*■) ~Z cly„(r) sin rn0/m 0 " r < m 

= 2A(/^ (clyQ(O) + 2^ ('Jy^(k) sin n rk/m) (1.191.b) 

S’xyg(r) is known as "the quadrature spectrum of y(t) and x(t)". Eqns, (1. 191. 
a-b) must be smoothed for the finite length to yield 5»y ( r) and 5*y ( r) . Finally 
we can write for the cross -covariance * ° 

|S;^(r)| = (Sx*y^(r)+Sxy^(r)l^ (1.192) 

and for the phase 

?^ ( r) = tan"^ ( S,y^( r)/Sx.v^ ( r) ) (1. 193) 

Now, we may summarize the main steps for the calculation of the cross 
spectrum: 

(1) Remove mean and trends. 

(2) Apply window function so that x(t) = x' (t)‘h( t) , y (t) = y*(t)*h(t) with 
the Fourier transforms X(f) » X'(f)* H(f), i'(f)=t'(f) 


AAA* 

<3) Compute tho Fourier ooefflolentt ei(f)i b>(f)( e.(f)t by(f) from ©qm, 
(1. 110) end (1. Ill), or better the ooefflolente (l«(k), 0 |(k), Ay(k), iiy(k)) 
from eiptB. (1.65.e-b) using FFT. 

(4) Compute Uie oospeotrum end qued-speotrum by reerrenglng eqn. (1. 183) Into 
reel end Imeglnery perts to get 

8«y,(k) - cii.(k) Sy(k)+bx(k) by(k)l/2 

A A A A A 

S*jj,(k) - t««(k) by(k) -ey(k) b*(k)l/2 

(5) Sum the resulte ebove over (2z 1) elementery bends to obteln 

r(Si<fl)4'i 

8»,(r| .(“■(k) a,(k)+b.(k) b,(k)|/2 

•“> _ '(C)- .... a-i”) 

S)iy-( r) “ I] f®»(k)by(k)-ey(k) bx(k)]/2 

0 kMr(*lTi)-t 

(6) Conqnite the cross spectrum lsl^( r) | and the fduse r) from eqns. 

(1. 192) and (1. 193) respectively. 


1,7, 6. 2 Coherence 
The followii^ expression 


1 + 

S*y(r) 

r-1 

.1- 

S‘xy( r) 

P 


(1. 196) 


is normally distributed with a variance approximately equal to (1/DF) (Rayner, 
1971, p. 98). Then the approximate limits of eqn. (1. 196) are given by 

arc tanh l|Sxy(r)ji ± (-”)^ (1.197) 

where (x = significance level, and = the expected limit, and v = degrees of 
freedom. 


1.8. Summary 

The computation of a spectrum from continuous and discrete data and 
spectral analysis has been explained in this chapter. In geodetic applications 
we generally have discrete and finite-length data. The selection of finite-length 
data causes undesirable effects on the spectrum known as the spectral 
leakage (Harris, 1978, pp. 51-52). In order to minimize the spectral leakage 
windows are applied as explained in this chapter. We believe that this chapter 
is a good reference for geodesists, who use discrete and finite-length data. 
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a. Spectral Analviii of QEOB»8 Altimeter Patti 
a»l Sattlllte Altlmetgy 


The determinetloD of the grevlty field end the shape of the earth has been 
the main goal of Geodesy for centuries, if we were able to measure gravity over 
the entire eartii, then a conslderaUe number of geodetic problems could be solved. 
However, In practice, we do not have enough observations In many parts of tiie 
earth, mainly In ocean areas. 

Recently the Geodynamics Experimental Ocean Satellite (GE06-3) obser- 
vations have made It possible to determine gravity anomalies from altimeter 
helf^t observations. By satellite altimetry we hope to measure the spectrum 
of the sea surface to a very high frequency [X (wavelength) 1000km ) so that 
we can determine small scale changes In addition to the Earth’s gravity field, 
lunl-solar effects, atmospheric tides etc. In fact, oceanographers would like 
to know sea topography at the 10 cm level from the mean sea level In order to 
detect ocean currents. 

The geometry of an altimeter borne satellite wltii respect to the oartii's 
center Is Illustrated In Figure 2. 1. Here we simply assume that the line OA 
connecting the geocenter and the satellite Is perpendicular to the Instantaneous 
sea surface (6S) and the mean sea surface (MSS). The altimeter measures the 
distance (h* *■ AP* ) between the satellite and the satellite's footprint of ISS. 

We can compute the coordinates of the satellite through the observations made 
at the tracking stations distributed in a world-wide network. Hence we can 
derive the geodetic (^,X) or Cartesian (X, Y, Z) coordinates of the satellite's 
footprint with respect to a reference ellipsoid as well as its distance (r' « OA ) 
from the center of this ellipsoid, and consequently the height (h AQ) of the 
altimeter above the ellipsoid. The sea surface height (SSH), which is the 
separation between the reference ellipsoid and ISS, Is computed as follows 
(Rapp, 1977, p. 2) 

SSH - h - (h' - a + b) (2.1) 

where R Is the refraction correction 

b Is the a priori altimeter bias (if any) 
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A <8ftt«mte) 



0 (Geooenter) 


Figure 2. 1. The Geometry of Satellite Altimetry. 


LH 


r»- 
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The geold la usually defined as the equipotential surface coinciding with 
MSS. The distance (N) between the reference ellipsoid and MSS is known as the 
geoid height or undulation. There exists a separation (^H) between ISS and MSS 
due mainly to tides, winds, storms, currents, etc. Therefore, SSH should be 
corrected for the factors mentioned above in order to obtain the geoid height (N) 
Hence we have 


N •= h-(h’ - R + b) -AH 


( 2 . 2 ) 


Eqn. (2. 2) is nothing else but an observation equation, the undulation N being 
considered as the observation. For one or more arcs Rapp (Ibid, p. 12) uses 
the following mathematical model for the adjustment of geoid heights: 


N+BA -N, + V 


(2.3) 


where N : 
1 » 
i * 

N* i 
Y 5 


column vector of geoid heights implied by altimeter date 

the design matrix of error model parameters 

column vector of parameters of the error model 

column vector of geoid heights implied by the reference surface 

(in our case GEM 9 surfoce to degree 20) 

column vector of residuals. 


For the adjustment purpose a data set of 419294 geoid heists implied by 
GEOS-3 altimeter data, in 2003 arc segments, was considered. First a primary 
adjustment was carried out, followed by regional adjustments. The details of the 
adjustment are described in Rai^ (Ibid, pp. 1-24). 
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2» 2 Owvity AnonuUy Raoowry 


As ws mentioned st the beflnoiof (A this ohspter, we sre mainly interested 
In the determination of point and mean gravity anomalies from the QE06-S altlm- 
elsr data. There are various teohnlques for the anomaly recovery. Hummel, 
SJoberg and Rapp (1277) have used <1) the inverse Stakes' equation or the 
Molodenskil equation (by Rummel); (2) the spectral analysis approach (by 
SJoberg) t (3) the least-squares collocattcn (by Rapp). In addition to these 
techniques, two new methc^s will bo introduced in this paper (1) the least- 
squares collocation using Toeplitz matrices; (2) the frequency domain least- 
squares collocation, in the third and fourth chapter respectively. 


In the next sections we will examine the power spectrum of the adjusted 
GE06-3 altimeter data and predicted point anomalies along various arc segments. 
Point anomalies have been computed from altimeter data by using least-squares 
collocation (Moritz, 1975) as follows: 

^ - C^.N (Cnn + P)~' N (2.4) 


where ^ ; 

n i 

t 


Snn S 

D 5 


predicted point anomalies 
vector of given geo id heights 

row-vector of the covariance between the anomaly predicted and 
the geoid heights 

square and symmetric covariance matrix of the given geoid heights 
the noise matrix of the given geoid heights (taken diagonal). 


The point anomalies predicted as described above do not reflect the contri- 
butions of the high-frequency components to the variance (total average power) 
due to the following facts; 

(1) Altimeter observations used in the adjustment procedure described 
above represent averages of two seconds in the low data rate, imposing a 28- 
km limit on the shortest wavelength of information. 

(2) We use a global covariance model (C^(n and Cnn above), which 
acts like a low-pass filter, in the least-squares collocation. 


The spectral analysis approach will enable us to point out the contribution 
of every existing frequency or wavelength component to the variance (total average 
power). Then we can make a reasonable conclusion about the resolution of altim- 
eter data and predicted anomalies. 


2. 3 Power Spectrum from Geoid Heights 

In this and following sections we will examine the breakdown of variance 
with respect to the frequencies (or wavelengths) and the contribution of each fre- 
quency conq>onent to the total average power. The power spectra of some ten 
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profU«t (numbcrtd oat to ton) tliown In Figure 2.2 Imve been oooiklered In thii 
Motion. From tho ndjustod 0E06-2 nltlmttor dnto wo hnvo oomputod goold 
helgbto along thoM proflloa at the Intoroootlooa of a latltudo-longltude giatloule 
of one dogoM or aio (Rappi 1970). Each profile haa been conaldered indepen- 
dently for the removal of the mean value along the particular profile from com- 
puted undulatlona and for the oomiMtotlan of the Fourier tranaform* which waa 
carried out by FFT algorlthma of IMSL llbrariee at 08U» by eqn. (1.02), l.e. 


X. ■ i"? X(k'4t) »•«. 


even 

odd 


(2.8) 


where X, ( the oonq>lex Fourier coefflc lento 

x(k*dt) t the geold height at the (k*dt)th point after the removal of the 
mean value along the profile 
N ; the mimber of undulatlona along the profile. 


In aection 1.6 we have defined the total average power (variance) of a 
fonction x(k*dt) ; k > 0, 1,..., (N-1) aa 


1 

P*v6* * Tf L U(k*At)l* 

*^11 «0 

in the time domain and by Paracval'a theorem an 
Pa.6* - E Pn* 


k MO 


( 2 . 6 ) 


(2.7) 


wtore 


Xo* for n 
2|xJ* for n 
for n 
zlXN^l^for 0 


0 

1,2,...,(Nh- 1) 
tin and N even 
Nm and N odd 


( 2 . 8 ) 


in the frequency domain. P^ , n - 0, 1, . . . , Nh la known aa the nth degree 
power of the periodic fonction x(k*At). Po ■■ Xo is the average value (alao 
called dc value) of Uie fonction x(k*At) . Since we have removed mean valuea 
before tranaformationa, they are equal to zero in all our computations. Even 
if there exists a mean value during the transformation, the power of zeroth 
degree la excluded from the variance to obtain 


Kvo* = PvG® - Po* = E Pn* 

n S2l 


(2.9) 


The contribution of the nth degree power to the variance is given by 


Cb - 



( 2 . 10 ) 


and the power implied by frequencies up to M < Nh (from now on we will call 
Mtk cumulative power) is given by 
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Figure 2.2. GEOS-3 GrouiK’ Tracks and Some Selected Sub- Profiles. 



(2.11.a) 


N 



finally, tbe contribution, of the M^*oumulatlve power la given by 

Rm « Sm/p 4 (2.11.b) 

We know that the frequency f« la defined aa 

f, - n/S (2. 12. a) 

where S la the length of the profile, and the wavelengHi la defined aa 

X. ■ S/n (2.12.b) 

Hence Pn^ can alao be called the power of nUi-wavelength component, and we 
can consider the wavelength aa the variable of the power. 

The ten profllea ahown In Figure 2. 1 have been transformed Into the fre- 
quency domain In order to obtain the degree power P„^, Its contribution c» to 
the variance, the M*'*-cumulatlve power to the variance. The profiles and the 
computation results are described In Table 2.1. 


Table 2. 1. The Breakdown of Variance According to Wavelengths. 


Profile 

No. 


Xj, g Xg 

Length 

Tl 3 



Wavelength 



Pav® 

X>6000 

X>3000 

O 

O 

X>1500 

X>1000 

X>500 

o 

o 

A 



_(Deg.) 

(km) 

_(m®) 


S 

cumulative power in m^ 



1 

-15,-15 

40,120 

8485 

342.8 

294.8 

318.1 

328.3 

331.9 

337.5 

339. 7 

342.8 

2 

-37,-37 

320,110 

13232 

367.6 

330.4 

357.1 

359.4 

361.1 

363.8 

365.9 

367.4 

3 

-48, -48 

10,208 

14658 

619.2 

577.9 

601.4 

608.4 

610.3 

614.0 

616.3 

618.6 

4 

-60, 60 

165, 165 

13232 

678.2 

603.4 

641.4 

652.3 

658.7 

665.8 

673.5 

678.2 

S 

-60, 65 

190, 190 

13899 

269.9 

160.8 

236.8 

245.6 

254.1 

258.0 

265.1 

269.4 

6 

-65, 65 

335,335 

14344 

521.8 

437.3 

487.5 

500.3 

504.7 

508.8 

516.2 

521.8 

7 

-60, 20 

60, 60 

8784 

1520.6 

1295.8 

1399.3 

1451.4 

1463.7 

1485. 6 

1505.4 

1520.6 

8 

20, 20 

120,250 

13479 

841.3 

756.2 

801.5 

816.8 

827.9 

830.1 

836.3 

841.3 

9 

-60, 65 

180, 180 

13899 

562.7 

458.9 

520.2 

534.2 

544.4 

550.9 

558.0 

562.7 

10 

Mean 

25, 25 

265,341 

7760 

713.2 

643.7 

605.8 

552.1 
1 

618.8 

598.2 

642.9 

614.0 

659.6 

621.1 

682.0 

629.7 

700.1 

637.7 

713.2 

643.7 


From the table above we see that the long-wavelength components of undula- 
tions contribute the greater part of the variance (power) P^, e. g. X>6000, 
X>3000, X>2000, X>1000, X>500, contribute 86 %, 93 %, 95 %, 98 %, 99 %, of 
the variance respectively. 
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2.4 Power Spectrum of Oeotd Hetehto with Respect to the GEM 9 Reference Surface 


Id the adjustmeat of GEOB-3 altimeter data a referoDoe aurface Implied by 
OEM 9 potential ooefficienta to degree and order 20 baa been iiaed. Therefore, 
we have decided to examine the powe r epectrum of geoid heii^ta with reapeot to 
the GEM 9 reference aurface, aay AN. M order to compute AN firat we have 
to compute reference aurface, i.e. GEM 9 surface, undulations and then subtract 
it from adjusted altimeter undulation at the same point. The GEM 9 surface un- 
dulations are computed as follows (Rapp, 1977): 


Nk 


GM 

yr 



cos mX <f ?£a ain mX) (sin o) 


(2. 13) 


where GM 


r 

a 



I the geocentric gravitati al constant 

; the geocentric distance . . the computation point 

: an equatorial radius 

; (iilly normalized potential coefficients 

: ftilly normalized Legendre ftmctions 

; the geocentric latitude and longitude 

: the normal gravity at <p . 


Thus we can write for the undulations with respect to the GEM 9 surface 


AN = N-N« (2.14) 

where N is the adjusted GEOS-3 geoid height at which Nk is being computed. 

Seven arcs of about 13000 km length have been selected for the spectral 
analysis. These arcs are described in Table 2. 2. 


Table 2.2. Arcs 2, 3, 4, 5, 6, 8, 9 and Statistics. 


Arc No. 

Latitude 

Longitude 

Length 

Mean (SU) 

RMS 

(See Fig. 2. 2) 

<Pi*<Pa (Deg.) 

» ^3 (DBR* ) 

(km) 

(m) 


2 

-37, -37 

320, 107 

13054 

0.3 

1.6 

3 

1 

00 

1 

00 

10, 186 

13021 

-0.6 

2.3 

4 

-60, 57 

165, 165 

13010 

-0.5 

2.1 

5 

-60, 57 

190, 190 

13010 

-0.6 

1.9 

6 

-65, 52 

335, 335 

13010 

0.7 

2.6 

8 

20, 20 

120,245 

13061 

0.9 

2.9 

9 

-60, 57 

180, 180 

13010 

-0.4 

2.6 
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2, 4, 1 Powr Speotnim when Data are Atwroed to be Pertodio 


tf we Meunie data are periodic, l.e. data repeat tiiemselvea Indefinitely, 
then we will have a dleorete apeotrum made up of a finite number of frequenclea. 
The geold helg^ta given with reapeot to the GEM 0 reference aurfoce, which are 
defined In eqn. (2. 14), have been tranaformed Into the frequency domain uaing 
FFT algorltiuna mentioned before. Finally the breakdown of variance, 1. e. 
the degree power, haa been computed by eqn. (2.8). The reaulta are described 
In Table 2.3. 


Table 2. 3. The Breakdown of Variance and Degree Powera with Reference to 


Wavenumber 

Correa- 

Wavelength 

Arc No. 

or 

ponding 

Global 

Degree 

(km) 

2 

3 


5 

6 

8 



9 

Mean | 

Frequency 

Degree Power 
(m ♦♦ 2) 

1 

( 3) 

13020 

0.17 

0.99 

0.20 

0.13 

0. 55 

0.18 

0.11 

0.33 

2 

( 6) 

6510 

0.44 

0.22 

0.03 

0.76 

0.44 

1.86 

0.54 

0.61 

3 

( 9) 

4340 

0.17 

0.30 

0.90 

0.54 

0.80 

1.31 

0.19 

0.60 

4 

(12) 

3255 

0.19 

0.52 

0.24 

..16 

1.46 

0.54 

0.95 

0.58 

5 

(15) 

2604 

0.33 

1.78 

0.51 

0.10 

0.41 

0.14 

0.15 

0.49 

6 

(18) 

2170 

0.08 

0.15 

0.05 

0.06 

0.03 

0.48 

0. 18 

0.15 

7 

(21) 

1860 

0.29 

0.01 

0.49 

0.00 

0.69 

0.01 

0.52 

0.29 

• 


0 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

• 

• 


0 

• 

1.67 

• 

3.97 

0 

2.42 

• 

1.75 

1 

• 

4.38 

• 

4.52 

• 

2.64 

c 

3.05 

a 

Pm 

i 

• 

• 

2.70 

0 

5.21 

• 

4.49 

. 1 

3.71' 

1 

1 

• 

6.76 

• 

8.48 

1 

• 

6.8l! 

0 

5.52 


From Table 2. 3 we see that the mean of variances of geoid heists Is 5. 52 
m^ . But there Is an uncertainty on each given geoid height, which is about ± 1 
m in our case. Then under the assumpticm of white noise, i.e. no correlation, 
we can compute the reduced power as follows (Wagner, 1977, p. 14) 

P ® (reduced) ■ P * (measured) - P* (noise) (2. 15) 

P® (reduced) = 5.52 - 1.00 = 4.42 


and finally P = 2. 1 m. 

So we have recovered about 2.1m undulation information in addition to the 
GEM 9 undulations. From Table 2. 3 we also see that X>2000 km contribute about 
1.6m , and X<2000km contribute about 1.3 m of this extra 2.1m undulation 
information. 

The results given in Table 2. 3 also yield a global power for harmonic degrees 















i ■ 3, 6, 9, . . . . However, we have to keep in mind dint we nnsumed periodicity 
of data and uaed only seven profiles, which are about 1/3 of a great circle around 
the earth. So in the next section we will consider non-periodic data and derive 
a global power spectrum. 


2. 4. 2 Power Spectrum from Non-Periodic Data 

In this section we are after a global degree power using GE06-3 altimeter 
data. However, we do not have any complete altimeter data around the earth. 

So we assume that (he power spectrum of each profile, which is a part of a great 
circle around the earth, is an estimate of the global power spectrum after scaling 
the frequencies and power. Since the variance of a sub-profile is a good estimate 
for the global variance, the degree power of each sub-profile contains about 
(Wagner 1977, pp. 89-93) 

Cl S (length of great circle)/ Si (length of sub -profile) (2.15) 
times the global degree power of the equivalent global frequency. 

Since the data of a great circle around the earth are periodic we can use 
a window to obtain a periodic function, say y ( n ), n =* 1 , 2, . . . , Nmx . from 
a non-periodic function, x(n), n = 1,2, . . , , N < Ni.x . We can explain the 
above procedure by considering the example above. We are given 120 observa- 
tions defined as 


x(k*Aj|2)), k = 0,1, ...,119 
and 

A0 = 1® 

Now let y( be defined as 

y(n«A((») = w(n) • x(n*<iV), (2.16) 

n = -120, -119, . . . , 0, . . . , 119, 120, . . . , 240 
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•uoh that 


^ ^ r 1 for 0 < n < 119 
" io otherwise 

then, y(n^0) Is defined along the whole great circle. So y(n*A0) is periodic 
and it can be transformed into the frequency domain and interpreted as global. 
However, in orter to maintain the variance (or average power) of y(n»L\l)) 
equal to the variance of x(k*A((>) , the degree powers of y (n should be 
multiplied by C| defined by eqn. (2. 15), which will be explained later. 

In order to compute a global power spectrum we have selected some four 
sUb-profiles with observations one degree apart, described in Table 2.4. 


Table 2.4. Arcs 4, 5, 6, 9 and Statistics. 


Arc No. 
(See Fig. 2.2) 

Latitude 
01 ,09 (Deg.) 

Longitude 

) 

Length 

(km) 

Mean (AN) 
(m) 

RMS 

4 

-60, 60 

165,166 

13343 

-0.5 

2.0 

5 

-60, 65 

190,190 

13899 

-0.4 

2.1 

6 

-65, 65 

335,335 

14455 

0.6 

2.5 

9 

-60, 65 

180,180 

13899 

-0.3 

2.6 


The global power spectrum of each sub-profile of Table 2. 4 has been com- 
puted as follows: 

(1) The mean of each arc was removed from the given data to obtain 

N-l 

x(k*A(^) = x(k*Zi(p) - x(k»Aj/)), k = 0, 1, . . . ,(N-1) 

where N is the iuLai number of observations along the particular arc and A0 » 1° 
in our examples. 


(2) x(k*A0) has been multiplied by a cosine -tape red window, which is 
explained in section 1. 7. 3. 1. 5, to obtain 



y (k*A((») = c (k*A0) • x(k*A0), 

n) s k n) 

(2.17) 

where 

1 

ill-cos (25TTk/N)l 

if 0 s ks n/25 



1 

c(k»At(') = j 

1 

if N/25 s k^24N/25 

(2. 18) 


i( 1-cos (25TT(N-k)/N)l 

if 26N/25<k «(N-1) 



0 

otherwise 
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and C <■ S (the length of great otrcle)/a (the length of eub-profile). 

(3) The variance of x(n*A(|))» n-0,l|..., (N~l) la defined by eqn. (2.6) 


aa 


Pavo* - — ^J^Jx(k'At)l* 


( 2 . 19 ) 


or equivalently 

P*vG* “ t* ly(h*^ll>)l*# where N,.x ■ C • N (2.20) 

Na«x kTO 

If we neglect the effect of the tiq>ering of the coalne-tapered window. The Fourier 
tranaform of y(k*A(()) haa been computed by eqn. (1.62), i.e. 

Yp « y y/lc»A((>) e"**’^'* (2,21) 

Naax kaO 

where Yp i the global complex Fourier coefficients (unacaled) 

Naax ; the dimenaion of periodic data vector y(k^V>) 

C ; the coefficient as defined above. 

Equation (2.21) is multiplied by in order to maintain the average power 
defined by eqn. (2. 19) or (2.20). 


(4) The breakdown of the variance according to wavenumber (or frequency) 
has been computed by eqn. (2. 8), namely 


where 


P« 


Y<j* for n “ 0 

2|Ynl for n‘^1,2,..., (Na~l) 

for n - N„ and Na.x even 
2 |Yn„|^ for n^N^and N..x odd 


Nmx/ 2 if Naax even 
(N..,-l)/2 if N„* odd 


( 2 . 22 ) 


and finally the variance (total average power) has been computed, 


Pavg* 


ns 0 


P 


2 

n 


(2. 23) 


The results are described in Table 2.5 for the sub-profile given in Table 2.4. 
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Table 2. 5. The Breakdown of Variance and Degree powere with Reference to 
OEM 9 Surfhoe from Non-Periodks Data. 


Wavenumber 

or 

Frequency 

Wavelength 

(km) 



Arc N< 

h 


4 

6 

6 


Mean 

Power 

(m**2) 

1 

40 

000 

0.01 

0.12 

0.15 

0. 12 

0.10 

2 

20 

000 

0.03 

0.21 

0.32 

0.19 

0.19 

3 

13 

333 

0.04 

0.12 

0.18 

0.04 

0.10 

4 

10 

000 

0.03 

0.16 

0.03 

0.06 

0.07 

5 

8 

000 

0.01 

0.39 

0.16 

0.41 

0.24 

6 

6 

667 

0.01 

0.34 

0.19 

0.67 

0.30 

7 

5 

714 

0.06 

0.05 

0.03 

0.62 

0.19 

8 

5 

000 

0.17 

0.06 

0.11 

0.48 

0.18 

9 

4 

444 

0.23 

0.27 

0.29 

0.41 

0.30 

10 

4 

000 

0.14 

0.19 

0.11 

0.35 

0.20 

11 

3 

636 

0.04 

0.01 

0.04 

0.24 

0.08 

12 

3 

333 

0.06 

0.09 

0.42 

0.13 

0.17 1 

13 

3 

077 

0.09 

0.20 

0.53 

0.06 

0.22 

14 

2 

857 

0.10 

0.12 

0.18 

0.05 

0.11 

15 

2 

667 

0.12 

0.04 

0.12 

0.06 

0.09 

16 1 

2 

500 

0.14 

0.02 

0.31 

0.10 

0.14 

17 

2 

353 

0.10 

0.00 

0.20 

0.16 

0.11 

18 

2 

222 

0.04 

0.04 

0.01 

0.16 

0.06 

19 

2 

105 

0.02 

0.11 

0.09 

0.09 

0.08 

20 

2 

000 

0.03 

0.06 

0.18 

0.00 

0.07 

21 

1 

905 

0.10 

0.01 

0.19 

0.05 

0.09 

22 

1 

818 

+0.14 

+0.03 

40.21 

+0.16 

+0.13 

• 

• 



1.70 

• 

2.64 

• 

4.05 

. 

4.61 

• 

3. 25 

• 

A *a 

Pa« 



• 

3.99 

• 

4.31 

• 

6.27 

• 

8.07 

• 

5.66 


As we see from Table 2. 5. the mean of variances (average powers) of 
undulations with reference to the GEM 9 surface is about 5. 66 m^ . If we con- 
sider aim uncertainty on each given data point, then with the assumption of white 
noise for the data, we can compute the reduced power as follows; 

"P® (reduced) = P ® (measured) - P®(noise) 

(reduced) = 5.66 - 1.00 = 4.66 m® (2.24) 

and fiirther P = 2. 16 m. 

Thus we have recovered about 2. 16 m undulation information in addition 
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to the OEM 0 unduUtionf* From Table 2*5 we alao aee that 2000km con* 
tribute about 1.7 m ^ and X < 2000km contribute about 1.3 m of thia extra 
2. 16 m undulation information. Notice the almoat identical reaulta given in 
Table 2. 3. where we had aaaumed periodicity, for reduced, meaaured powers 
and extra undulation information. 


2, 5 Power Spectrum of Gravity Anomalies 

As we mentioned before the main goal of Geodesy is the determination of 
the gravity field and the shape of the earth. In sections 2, 1 and 2. 2 we have 
e:q)lained how to detain geoid heights from satellite altimeter data and how to 
recover gravity anomalies from these geoid heights, hi this section we will 
examine the power spectrum of these recovered anomalies. We are mainly 
interested in pointing out the resolution of these anomalies. In other words, 
we will try to point out Uie minimum anomaly wavelength recoverable from 
these anomalies. 

In our computations throughout this section the N-data vector has been 
transformed into the frequency domain by eqn. (1.82) after the removal of mean 
and discontinuities at tiie ends. Discontinuities at the ends have been removed 
by a cosine-tapered window. Degree powers P„^ and the total average power 
have been computed from eqns. (1.80) and (1.81) respectively. Frequencies 
and wavelengths have been computed from eqns. (2. 12. a) and (2. 12. b). 


2. 5. 1 Power Spectrum from Recovered Anomalies Only 

Two GEOS-3 arcs, one in the Calibration area and the other in the Philip- 
pines area, have been selected. These arcs are shown in Figure 2. 2. Along 
these arcs point anomalies have been computed from adjusted undulations at 
the data points using least-squares collocation. The data vectors of point anom- 
alies have been transformed into the frequency domain in order to obtain degree 
powers, their contribution to the variance, and finally the cumulative contribution 
by eqns. (2.8), (2.10) and (2. ll.b) respectively. The results are described in 
Table 2. 6. 

The altimeter data used in the adjustment are approximately 2-second 
averages. So the shortest wavelength of information recoverable from the ad- 
justed GEOS-3 altimeter data of OSU is roughly 28 km if we assume the data are 
noise-free, fa fact we know that they are not noise free. These undulations 
have about 1 m standard deviations. 
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Table 2, 6* The Breakdown of the Variance of Recovered Point Gravity A nom alies. 
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From the remits given In Table 2.6 we can point out the shortest anomaly 
wavelength recoverable from point anomal lea implied by GEOB>3 altimeter data. 
First we will asmme ttie standard deviations of point anomalies to be about 27 
mgals, this having been found In the solution of least-squares collocation. Second- 
ly we will also assume that the noise spectrum is equally distributed on each degree. 
Then we have, with (1.07) 


, n % 172® nn^l*/268 = 2.72 mgal^for the U'*' arc 

^ " " 1729 mgnlVsOO » 2.43 mgal* for the 12“’ arc 


(2.25) 


We see from Table 2.6 that tiie si^ial -to -noise ratio is greater than one 
above the 264 km wavelength for the 11^ arc and above the 115 km wavelength 
for the 12*" arc. 


2. 5.2 Power Spectrum from Measured (by Ship) and Recovered Point Anomaliea 

In the previous section we have examined the power spectrum of computed 
point anomalies along two arcs separately. In this section we will examine the 
power spectrum of anomalies using predicted point anomalies and measured point 
anomalies, which were kindly provided to us by Mr. W. H. Chapman, U.S. Geological 
Survey. The anomalies are along the sub-profiles 13, 14, and 15, which are de- 
scribed in Table 2.7 and shown in Figure 2.3. 


Table 2.7. Statistics of Recovered (Predicted) and Measured Point Gravity 
Anomalies Along the Arcs 13, 14, 15. 


m 

N 

Length 

Latitude 

Longitude 

Mean 

Variance 

Mean Diff. 

RMS 




Xi iXa 

(mgals) 

uSIm! 

**2) 

(mgals) 

Diff. 

■ 


(km) 

(Deg.) 

(Deg.) 

m 

Meas. 

Pred, 

Meas. 

(Meas-Pred) 

(mgals) 

13 

132 

383 

25.06, 22.36 

122.11, 124.45 

-43.8 

-45.4 

3782. 

6776. 

-1.6 

39.9 

14 

78 

203 

22.54, 24.39 

124.53, 124.71 

-58.5 

-49.2 

4529. 

4109. 

9.3 

23.4 

15 

168 

473 

18.08, 13.88 

127.15, 126.28 

17.1 

13.9 

372. 


-3.2 

13.0 


The data vectors of predicted point anomalies and measured point anomalies 
have been transformed into the frequency domain separately along each sub-profile, 
given in Table 2.7, in order to compute degree powers, their contribution to the 
variance, and finally the cunmlative contribution by eqns. (2. 8), (2. 10), and (2. 11. b) 
respectively. The results are described in Table 2.8. 
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Tiible 2. 8. The Breakdown of Variance and Degree Powers of Point Anonalies Along Arcs 13, ii, 15. 
arc NC : 13 ARC NO : 14 ARC NO : 15 
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If we ftMume the eteiKkrd devletlonc of predicted anomellee to be 27 mgali 
ee la aeotloo 2. 5. 1| then we oaa calculate the uacertalntiea on each degree fiowcr 
to be 

729/132 ■ 5«52ntgal* for the 13*'' arc 
var(PRVri« 78 *9.35mgal” for the 14'* arc <2.26) 

729/148 ■ 4.92mgal* for the IS** arc 


We aee fitmt Table 2,8 that for predicted anomallea the signal >to-nolae 
ratio la greater Uian one above the 127km wavelength for the 13** arc and above 
the 68 km wavelength for the 14** arc, and finally above the 208km wavelength 
for the 15*" arc. The 208 km shortest recoverable wavelength of the 16'* arc is 
much higher than 127 and 68 km of the 13** and 14** arcs respectively. Assuming 
the same precision for predicted anomalies along each arc we can state that the 
big difference above Is most likely due to Luiooth ocean bottom topography along 
the IS** arc and rough ocean bottom topographies al<mg the 13** and 14** arcs. 


Thus, after all these analyses, we can conclude that the shortest point 
anomaly wavelength implied by GE06-3 altimeter data is about 70 km when the 
standard deviation of predicted anomalies is about 27mgals. 

Now consider the standard deviation of each measured point anomaly to be 
Smgals. Then we can compute uncertainties on each degree power as follows: 

25/132 = 0. 19 mgal® for the 13** arc 
26/ 78 = 0. 32 mgal® for the 14'* arc (2. 27) 

25/148 - 0. 17 m^l^ for the 15** arc 

Ue sec from Table 2. 8 that for measured point anomalies, the signal-to- 
noise ratio is greater than we above the 12. 3km wavelength for the 13** arc 
and above the 6.3 km wavelength for the 14** arc and finally above the 21.5 km 
wavelength for the 15** arc. We conclude from the results above that the meas- 
urements (ship) have a much greater high-frequency contribution to the power 
than predicted point anomalies; this we had anticipated beforehand. 


var(PR)M«>, =» 


2.5.3 Power Spectrum from Measured (by Ship) and Recovered Point Anomalies, 
Bathymetry, and Filtering 

m previous sections we have examined the power spectrum from predicted 
anomalies, as well as measured (by ship) and predicted anomalies. In this section 
we have selected three more profiles, numbered 16, 17, and 18, which are shown 
in Figure 2.3. Along these, measured point anomalies and bathymetry were 
provided to us by A. B. Watts of Columbia University. On these profiles we have 
computed point anomalies using least-squares collocation from CEOS-3 altimeter 
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Figure 2. 


, GEOS-3 Ground Tracks and Some Selected Sub-Profiles. 
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data. Predlctiona have been made at those points where measurt>d anomalies have 
their local maxima or minima and are still not hirther than .*10 km from m^lKhlxiring 
prediction points. These three profiles are described In Table 2.9. 


Table 2.9. Statistics of Recovered (PrecllekHl) and Measured Point (iruvitv 
Anomalies Along the Arcs 16, 17, 18, 


aa 

HP 


Length 

Latitude 

Longitude 

Mean 

Variance 

Mean 


HMS 



(km) 

Oi ,«a 

X\ fXa 

(mgals) 

(mgal^) 

Diff. 


Diff. 

■ 

m 



(Deg.) 

(Deg*) 

Pred. 

Meas. 

Pred. 

Meas. 

(Mcas-Pred) (mgals) 

16 

62 

390 


- 3.19,- 4.00 

85.12,93.52 

-19.5 

-18.3 

689. 


- 1.2 

11.5 

17 

54 

178 

msi 

-16.76,-17.51 


-24.5 

-31.0 

1240. 

1791. 

6.5 


10 . 7 

18 

66 

210 

1013. 

-20.53,-20.64 


-8.1 

- 22.2 

153. 


14.1 


14.7 


We have also computed geold undulations by least-squares collocation at 
those points where point anomalies are predlckd along the profiles 16, 17, and 
18, All these, the bathymetry, measured anomaly, predicted anomaly, and 
geoid undulation profiles (corresponding to the profiles 10, 17, and 18) are plotted 
in Figures 2.4, 2.5, and 2.6 respectively. 

The data vectors of predicted |X)lnt anomalies and measured point anomalies 
have been transformed Into the frequency domain along each sub-profile, given in 
Table 2.9, in order to obtain degree powers, their contribution to the variance, 
and finally the cumulative contribution by eqns. (2, 8), (2. 10), and (2. ll.b) 
respectively. The results are described in Table 2. 10. 
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Table 2.10. The Breakdown of Variance and Degree Powers of Point Anomalies Along Arcs 16, 17, 18, 
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Af we did In previous seotlous let us sssume that the standard deviations of 
predicted anomalies are 27 mgals each. Then we can compute uncertainties on 
each degree powert 


! 729/62 ■ 11.76 mgal* for die IB^^arc 
729/54" 13. 60 mgal” for the 17^'*arc (2» 28) 

729/66" 11.05 mgal” for the 18^''arc 


We see from Table 2. 10 that for piredloted anomalies the slgnal-to-nolse ratio 
Is greater than one above the 200. 3km wavelength for the 16^** arc and above the 
125.9 km wavelength for the 17^** arc and finally above the 251.0 km wavelength 
for the 18^ arc. 

Let us again consider, as we did in previous sections, the standard deviation 
of each measured point anomaly to be roughly 5 mgals. Then we can compute 
uncertainties on each degree power as follows; 

{ 25/390 = 0.06 mgal” for the IB^^arc 
25/178 = 0. 14 mgal” for the 17‘^arc (2. 29) 

25/210 » 0. 12 mgal” for the 18*"arc 


We see from Table 2.10 that for measured anomalies, the signal-to-noise 
ratio is greater than one above the 22. 3 km wavelength for the 16*" arc and 
above the 11.4 km wavelength for the 17*" arc and finally above the 17.6km 
wavelength for the 18*" arc. 

We see again that measurements reflect greater high-frequency contribution 
to the total power than predicted anomalies. 

Since measurements reflect the shorter wavelength features of anomalies, 
the RMS difference between predictions and observations along a particular profile 
should get smaller when we use a low-pass filter for: 

( 1 ) measurements 

(2) both measurements and predictions. 

The low-pass filter used here can be described as follows: The data vector 
of the Mp 'tial (or time) domain representation along a particular arc has been 
transfc rmed into the frequency domain. Then the Fourier coefficients have been 
truncated above M ^ Nh , where Nh is the Nyquist frequency and M is the cut- 
off frequency in tiie basic interval, i.e. the length of the arc in question. Trun- 
cating the frequency domain series above the frequency M corresponds to low- 
pass filtering. After having generated the complex Fourier coefficients Xn, 
applying a filter only involves the truncation of the coefficients above M ^ Nh . 
Then by an inverse Fourier transform we can regenerate the spatial domain 


- 66 - 


repveieatation. Ja order to make our approach more underatandable let ua give 
an example: Suppoae we have an arc of 1000 km length with 100 equally spaced 
obse^rvationa ( N ■■ 100) on it. Then the Nyquist frequency la: 

Nh - N/2 - 50 (in the interval of 1000km) 

and the frequency, with eqn, (2. 12. a), la 

f« ■* n (in the interval 1000 km), n»0,l,...,S0 

and finally the wavelength, with eqn. (2.12.b), is 

X n ** l/fj^ * 1000 km/n , n * 1, 2, . . . , 50 


Thus we can compute the shortest wavelength recoverable, which corres- 
pmuls to the Nyquist frequency, for this particular example, 

Xn^“ 1000 km/60 “ 20 km 

The data vector of observations, say x ( t ) , t » 0, 1, ...» 99, is transformed into 
the frequency domain by eqn. (2. 5), 1. e. 

•t 

Xn * x(k)e’‘*”‘‘"/^ n =0,1,.. .,50 

The inverse Fourier transform (c.f. eqn. (1.71) ) yields 

x(k) = Xi(k) + Xi*(k) , k= 0,1,. ..,99 

•e 

where Xi(k) = Z„ such that 

Am = 50 km = 1000 km/M M = 20 

Now suppose we want to filter the powe r c(mtribution of any wavelength 
smaller than 50 km. Then we have 

X,^ = 50 km = 1000 km/M -* M = 20 

Now the truncated complex Fourier coefficients, say Y#, are defined as 

Y = n =0,1,..., 19, 20 

" ^0 otherwise 

Then, finally, we can regenerate the spatial domain representation, say y ( t ) , 
corresponding to the frequency domain representation Yn . The function y ( t ) 
does not have any high-frequency components above 20. The definition of y( t ) 
is as follows : 


Z ^ JX„ for OS n s 50 
" “ \ 0 otherwise 
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yi(t)+yi*(t) 


where 


y(t) ■ 

tf 

yUM ■ E 

^ n M n 


By ueing low -pan filtering aa deaoribed above we have filtered out any 
power contributiona from wavelengtha below 60 km, 100 km, and 200 km for the 
data veotora of meaaured anomaliea and bathymetry along the aub-arca 16, 17 
and 18 ahown in Figure 2.3, The regenerated apatial domain repreaentatione 
are ahown in Figurea 2.7, 2.8, 2.9 for X & 50 km, X ^ 100km, X ^ 200 km, 
reapeotively. In addition to the three eub-profilea above we applied low-paae 
filtering to the data vector of meaaured anomaliea along the eub-profiles 13, 

14, 15 of Figure 2.3 for X < 50 km, X £ 100km, and X ^ 200km. The 
reaulta are deaoribed in Table 2. 11. 

In the procedurea above we have only filtered out the high-frequency com- 
ponenta of the power of meaaured anomaliea. In addition we aiao filtered out the 
high-frequency componenta of the predicted anomaliea. The previoualy defined 
low-paaa filter haa been applied for the data vector of meaaured anomalies, 
bathymetry, aa well as predicted anomaliea along the profiles 16, 17, and 18. 

The regenerated spatial domain representations are shown in Figure 2.10, 

2.11, and 2.12 for X ^ 50 km, X ^ 100km, X ^ 200 km respectively. We 
also applied the low-pass filter to the data vector of measured anomalies and 
predicted anomalies along the sub-profiles 13, 14, 15, of Figure 2.3 for X^SOkm, 
X ^ 100km, and X ^ 200 km, The results are described in Table 2.11. 

Table 2.11. RMS Differences when Low-Pass Filtering is Applied to Point 


Anomalies Along the Arcs 13, 14, 15, 16, 17, 18. 


Arc 

No. 

Mean Diff. 
(Pred-Meas) 
(mgals) 

RMS 

Diff. 

(mgals) 

RMS when Xs 50 km 
filtered 

RMS when X^ 100 km 
filtered 

RMS u 
f 

'hen X^ 200 km 
lltered 

Meas. 

Both 

Pred&Meas. 

Meas. 

Both 

Pred& Meas. 

Meas. 

Both 

Pred& Meas. 

13 

6.5 

41.7 

40.4 

40.3 

32.0 

31.5 

18.9 

18.9 

14 

0.0 

26.6 

23.5 

23.2 

17.7 

17.4 

17.5 

12.7 

15 

0.9 

13.6 

12.6 

12.7 

11.8 

11.7 

8.6 

8.4 

16 

-1.2 

11.5 

11.7 

11.3 

11.9 

10.9 

10.3 

8.2 

17 

6.5 

18.6 

15.3 

15.1 

13.0 

12.5 

13.2 

12.0 

18 

14.1 

14.4 

13.1 

13.1 

11.1 

11.1 

8.6 

8.2 

Mean 

4.5 

21.1 

19.4 

19.3 

16.3 

15.9 

12.9 

11.4 


We see from Table 2. 11 that the RMS differences between the cases ( i) 
when only the data vector of measured anomalies are filtered and (ii) when 
both data vectors of measured anomalies and predicted anomalies are filtered. 
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Figure 2. 7. The Low-Pass Filtering 
(XsSOkm) on Measured Anomalies 
and Bathymetry Along Arcs 16, 17, 18. 

1. Altimeter Anomaly (without filtering) 

2. Measured Anomaly (low-pass filtering) 

3. Bathymetry (low-pass filtering) 
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Figure 2. 11. The Low>paee Filtering 
(X< 100km) on Altimeter Anomaliee 
Meaaured Anomaliee and Bathymetry 
Along Area 16» 17, 18. 

1. Altimeter Anomaly 

2. Meaaured Anomaly 

3. Bathymetry 







Are veiy eniAll* The dIflerenoeA are AinAll bectuie the datA vector of predlctlone 
refloota very llttiie power contributlona of high frequenoiei . We alao aee from 
Table 2. 11 that the mean of RM8 dlffereocea between predicted and meaauK^d 
anomallea la 21.1 mgala out of 6 aub*profUea 13, 14, 15, 16, 17, and 18. it la 
about 11.4 mgilB when the power contributlona from wavelengtha telow 200 km 
are filtered out. 
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8, Simple and Block Toaplita Mfctrioei 


8» 1 Simple Towlitg Matrioai 


A ToeplltB matrix la an (N***!* N^-l) matr'x ■ tkj • where tio > tt,.j 
k, J ■ 0, 1, . . . , N» that la to aay that t^ j la a ftmctlon of ( k-J ) rather than of 
(k,)) aeparately. So Tn la a matrix of the allowing foi*m 


Tn 


to t»i t«>9 • • • t~K 

tx to t— j * • • t»/N»t\ 

• • • • ' ' 

• • • • 

• t • e 

e ^ e e 

l»i til-xlN-* to 


(3.1) 


Covariance matrices of weakly stationary stochastic time series, and matrix 
representations of linear tlme-lnvarUint discrete time filters are of Toeplitz 
form (Gray, 1977). 


The very special structure of the Tn matrix can be e;^loited during the 
process of inversion, which yields significant savings in computational time and 
storage, to the past, several authors [ (Levinson, 1947), (Trench, 1964), (Kutikov, 
1967), (Zohar, 1969), etc.] have approached this problem and given algorithms to 
carry out the inversion of general Toeplitz matrices. 


Since we generally deal with ’positive definite and symmetric normal equations' 
in least-squares adjustment and/or collocation, only the inversion of symmetric and 
positive definite Toeplitz matrices will be considered here. 


Let us now consider the least-squares collocation model 


X = AX + s'+n (3. 2. a) 

which yields the following e;q)ression for the s-signal vector (Moritz, 1975, p. 15) 

S = Cf xS (X— AX), C ~ Cf* ,• + Cnn (3. 2.b) 

or equivalently, when the trend is removed before hand 

s = C.X X (3.2.C) 

where s is the p-signal vector, x is the (N+l)-measurement vector, X is the 
u-parameter vector with the design matrix A , C,x is the [p, (N+1) ] -crosscovari- 
ance matrix between 8 and X, Cx'a* is the covariance matrix between s' signals, 
Cnn is the covariance matrix of the noise, and C is the [(N+1), (N+l)j covariance 
matrix of the measurements x. If the number of observations is large, then tlie 
classical inversion of C in eqn. (3. 2) causes the greatest problem. The 5 matrix 
is of Toeplitz form if we assume* 
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> obtexvfttiont are at equally apaoed pointa along a profile 

- nolaea are equal 

- the covariance function uaed in Uie computation of C haa the property 

)jk ■ (C^t* )j-k * ^1* impliea the atationarity of the covariance 
Amotion. 

Witii the aaaumptiona above, it followa that: 

(1) Tn ■ S la non*aingular 

(2) Tn ia aymmetric, i.e. tkj ■ tjk 

^ 3 ) ^ Hk tkj aj ^0, n**0,l,...,fi for any aet a^ , m n 

" such that at leaat one of the ia not zero and finite. 

A reoumive procedure for the Inveraion of symmetric and positive definite 
Tn will be given here. This recursive inversion requires numerical computations 
proportional to (N-tl)* compared to (N-tl)^ of classical bordering methods and 
the storing of only one row ot the input and the output 


3.2 inversion of Simple Toeplltz Matrices 


It is well known that the inveroe of a symmetric and positive definite matrix 
is also symmetric and positive definite. Then it is possible to express the inverse 
Tn~^ of Tn as (Faddeyev, 1963) 


Tn’^ = a B 


where 


1 0 


B * 


bio I 
bao bgj 


0 

0 

1 



Ono bna* ..... 1 



(3. 3) 


(3.4) 


In complete agreement with the equality above we can write (Kutikov. 1966) 


A*^ = B Tn B^ 


(3.5) 


or equivalently, 


n 

X J “ ® » 

kiS^O 

b„„ - 1, 


^ bjk * tkn 


a, 


.-a 


na • 


kmo 


If we substtiute the following condition 


J = 0,1,..., (n-1), 
n “ 0, 1, . . . , N 
n» 0,1,. ..,N 


n » 1,...,N 


(3. 6. a) 
(3.6.b) 
(3.6.C) 


tkj * tic— J * h, j = 0,l, ..., N 


in the above equations, we obtain 
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m 

“Of J * 0| If • • • I n*lf«««fN 

If n • Of If • • • f N 

■ <if» n-0,l,...,N 


k «e 

b., 

a 


kaO 


(3. 7. a) 
(3.7.b) 
(3.8) 


Since Tn is non-degene mte by definition, b.k values are uniquely defined by 
the relations above. We will show that (Ibid, p. 68) under conditions (1), (2), (3) 
we have the following recurrence relatimist 


bn a “ 

t, 

n 

■ Of If.i 

,..N 

(3. 9. a) 


* bnii*t|j.f^f 

n 

* Of 1;> • • 

..(N-1) 

(3.9.b) 


^na “ bni(*tn^ 

n 

* Of 1 f • • 

..N 

(3.9. c) 

bftf 1^0 

■ “©n* ^nn • 

n 

® Of If t • 

..(N-l) 

(3.9.d) 

bn^lfk 


n 

* l#2f .. 

.,(N-1), k»l,,, 

.,(n-l) (3.9.e) 


In Older to show this relation let us write the left hand sides of eqn. (ii, 7,a) replacing 
n by (n+1) : 

n + 1 

• 8“0,l,...,n (3.10) 

and let us also assume that for the b,k k,m « 0 n defined by conditions (3.9) 

the relations (3. 7)-(3. 8) are satisfied. Using eqn. (3.9) we can writo 

r, 

X| “ b„+i,o*t-» + bn+i,k*t|{_, + bn+i,„4.i*t„+i.i 

~ bn+l,0*t_, + y bn,k-l *tk«» + ba+i^O V ba,a_),*tk.i + b„+ija4.1*ta..l_, 

•‘'*1 . k^l 

® [t_, ^ b 

k^x ■ — . « 

uBliag the equality (3. 9. a) and substituting m = k + 1, we obtain 

n-l 


n9n-k*t)(^] + [ \ bn,*-i*t,.f + • t„^l« g ] 

■^1 


“7 . 

X * bn+1,0 [bo,ot-. + y ba,„.k*tk-,l + [V b„,^ tk-,+ 1 +ba,n tn+i-,] 
jL, 1 » kS’o 

i-K *tk~» + y bn,k *tk-»+l f S=0.1,.,,,n 

k = 0 


X = b 


(3.11) 


By virtue of eqn. (3. 7, a), the right-hand sides of the above equalities are equal 
to zero for s = 1 n and by (3.9.b)-(3.9.d) for s = 0 also. 

Thus when we consider eqn. (3. 7. a), then the eqns. (3. 7) -(3. 8) follow from 
the conditions (2. 9) by induction. Since b„k n *= 0, 1, . . . , N a.id k = 0, . . . , n are 
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unique, we conclude thnteqne. (3.7) and (3.8) are equivalent to eqn. (3.9). 


liet k,J " 0, 1, . . . , N denote the elements of the Inverse Ti7^ such that 


Tn 


.1 


roo ro“i 

r{L rfi 




n»x 


Ton 

r^N 




(3. 12) 


then the following equalities hold (Kutikov, p. 66) 

tkj • fill + h„., 1,1, ^♦1/1+1 1 ^ k, J “ 0, 1, . . . ,(N+1) (3.13) 

This follows from eqn. (3.3), where Is taken to be equal to zero. 


Since the Toeplitz matrix Tn satisfying the conditions ( 1 ), ( 2 ) and ( 3 ) is 
symmetric and persymmetric (i.e. it has symmetry both its diagonal and cross- 
diagonal), in other words 

^kj ~ tji, , t|,j “ t<-j,N~k 

we can show that the inverse Tn'^ is also symmetric and persymmetric (Zohar, 
1969). The symmetricity of the inverse is apparent, so, we will only show that 
Tn"^ is persymmetric. For the proof let us introduce the Exchange matrix E 
defined as a square matrix with units along the cross -diagonal, and zeros else- 
where. Now suppose A is an arbitrary (NxN) matrix. If we examine the 
matrix product ( E E ), we see that the overall effect Is to exchange elements 
of A which are located symmetrically with repsect to the cross-diagonal. So 
a persymmetric matrix is a matrix H satisfying E H ^ E = H . Now from this 
formulation it is easy to show that TiT^ is persymmetric. Let us start with 

(Tn-^) (Tn') = I 
and note that E E - 1 to obtain 

E(tJ)"NeE)(Tn ) E = EE = I 
(E(tJ)"' e] (et; E] = I 


but ETnE = Tn. Hence IE(Tn)"^EJ “ Tn"’^ and Tn"^ is persymmetric. 

Toeplitz inverses have a very useful property; If the last row of the inverse, 
say rN*j, j = 0, 1, . . . , N , is known, then all other elements of the inverse can be 
computed iteratively from rN*!) . First let us compute rN^ by writing from eqn. 

N (N-1) . a® K 

Tkj “ +bN k «><,j 
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D , tvi > 1 , we obtain 


and letting k*N and alao noting 

ri?,j ■ am b«,j J "Otlt.t.tN 


(3.14) 


In order to compute the remaining elementa of the Inverae we have from eqn. 

(3. 13) 

(3.16) 


N N-1 8 

*NM W,J«1 


and alao conalderlng peraymmetry 

N N N-i a . 

TicJ * ni|-l,N.k “ -k buyii.j 0NH 

N **X 3 

“ am bN^_k (3.16) 

N * I 

Substitute the equality of from eqn. (3. 16) In eqn. (3. 15) to obtain 

N ^ u a u u a 

“ Tkj " dn.n-j am b>i|,N_i, + bN.n_x a^N (3>17) 


By conalderlng symmetry and peraymmetry, the algorithm can be summarized 
as follows: 


(1) Compute the last row of the B matrix and through the recursive 
procedure defined by eqn. (3.9); 


bnn " 1» 
«n 




* fk+x» 

t bnk * tn-k » 

0 

bn+X,0 " “®n*a|ju, 

bn^lfj * ^ bu+x.o b^, 

'^w-l,n-J+l “ bn^n-J bn+i^obn,k-i» 


n * 0,1,. ...N 
n “ 0,1,..., (N-1) 
n ^ 0, 1, . . . , N 
n = 0,1,..., (N-1) 

n = 1,2,..., (N-1) 

J - l,2,...,((n-l)/2) 


(3.18) 


By the formulae above the reduction In computer storage Is considerable, because 
bn+i,j, bn+i,n-i+i following the calculation can be stored In the place of bB,j«i , 
bfi^n -k • which are not used in any further calculation. 

(2) Compute j = 0, 1, . . . , N and the remaining elements of the Inverse 

byeqns.(3.14) and (3.17); 
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N 


N M N N 

rjN ■ ■ foji-i • n«} from •ymmetry and 

pertymmetry 


N (3* 


19) 


■ *\j “ b>i,N-ii**5i*bi(,N-j + ♦ k ■ (H4 + 1)» * • ♦ tN» J ■ (N-k), , . . ,k 

mihmrm M . /(N^J'/a if N odd 
whtro Nm - I' I, 


finally “ tS-k+i,N-j*t - from aymmetiy and per- 

aymmetry. 


Aa we nee from the* eiq[>re8aiono above, we do not have to store all the 
elements of the inverse at once. As soon as an element of the inverse is com* 
puted it is multiplied by ihe corresponding element in x-vector of our model 
(3. 2.2). So we only have to store one row of the inverse in order to solve the 
recursive equation <3. 17). A computer program written in FORTRAN language 
is given in Appendix 3. A. 


Using this algorithm some test runs have been made in double precision 
arithmetic on Amdahl-470 computer. The results are given in Table 3. l.for the 
solution of the linear equality 

Tn F - X (3. 20) 

where T,, is a [(N't-!), (N-i-l)l Toeplitz matrix whose inverse has been computed 
row by 7.*ow, X is the data vector of length (N 1) and F is the solution vector 
of length (N * 1) we dee ire. 

Ttd}le 3. 1. CPU Time and Core Storage for the Inversion 
of Simple Toeplitz Matrices. 


Dimension 

(N+1) 

CPU Time 
(sec.) 

Storage 

(K) 

100 

0.12 

3.2 

200 

0.48 

6.4 

400 

2.86 

12.8 

1000 

18.21 

32.0 


From the Table 3. 1 we see that the simple Toeplitz inversion is very efficient 
as far as time and storage are concerned. If we have equally spaced data with 
equal noises, then the covariance matrix of the data is of Toeplitz form and we 
obtain identical solutions with the classical inversion algorithms. But if the ob- 
servations are not equally spaced and do not have equal noises, then we have to 
use interpolated data and give equal noises to the observations in order to use 
the algorithm above. Thus we introduce some approximations to our solution. 

We can proceed and use this algorithm if we can show that the degree of approx- 
imation has negligible effect on the solution vector. To demonstrate the efficiency 
of this algorithm and to see the effect of approximation on the solution, two arcs 
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0 


of altinittry deiorlbad In Table S.2 and ahoarn In Figure 2.2 from QE08*3 have 
been examined. 


Table S. 2. The Deaorlption of Aroe 11 and 12. 


Arc 

N6. 

Area 

Latitudes (Deg.) 

Longitudes (Deg. ) 
Xi Xa 

Length 

(km) 

Nb. of Observations 

U 

12 

Calibration 

Philippines 

12.8480 

-0.0241 

30.0526 

32,4656 

300.0578 

146.5301 

201.0061 

126.8220 

H 

266 

300 


Aa the firatatep, the point gxnvity anomaliea hat'e been predicted at obaervation 
pointa by equ. (3. 2.o) uaing the original obaervationa a^ their uncertaintiea. Aa 
the aeoond atep, equally apaoed data have been created from the original obaer- 
vationa and the mean uncertainty of the original obaervationa along Uie particular 
arc haa been given to each predicted data point aa the uncertainty. Aa the third 
atep the point gravity anomaliea have been computed at the original obaervation 
pointa by eqn. (3. 2.c), but uaing the aimple Toeplitz Inveraion algorithm to 
invert C . Hie atatiatlca of thin atudy are given in Table 3. 3 and anomaliea 
are plotted in Figurea 3. 1 and 3. 2 for the 11^ and 12%rc reapectlvely. 


Table 3. 3, The Statlaticn of Aren 11 and 12. 


Arc No. 

RMS 

Max. Diff. 

Mean Diff. 


(mgals) 

(mgals) 

(mgals) 

11 

1.1 

5.5 

-0.1 

12 

4.8 

32.4 

-0.3 


The average atandard deviation of the predicted point anomaliea in about 
27 mgala. The RMS differences of 1. 1 mgals for the eleventh arc, and 4. 8 mgals 
for the twelfth arc above are much small' than the standard deviations of the 
predictions. So it would not be unfair to say that the approximations introduced 
in the examples above have negligible effects on the computed point anomalies. 

The twelfth arc of the Table 3. 2 was also used for ( l°x 1°) mean anomaly 
recovery of some 25 blocks, which are located in a ( 5*^x S") block whose co- 
ordinates are: 

(Pn= 10®, tPs = 6®, Xe = 137®, Xw = 132° 

The RMS difference between a rigorous least-squares collocation solution and 
Toeplitz solution has been found to be only 0. 1 mgals with a maximum difference 
0.2 mgals . The difference between the two solutions is negligible should we 
consider the standard deviations to the predictions, which is about IS mgals. 
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Fi^re 3. 1. Altlmetar Anomftly Using Toeplits Algorithm and Conventional 
Collocation Along Are 11* 
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Figure 3.2. Altimeter Anomaly Using Toeplitz Algorithm and ^ Toeplitz Algorithm 
Conventional Collocation Along Arc 12. 
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8.8 Blook-TonllU Itetricw 

Wi hAv« <Mflntd 1ti« •Imple To«pUts matrix m a matrix Tn wboaa (kj )di 
atommit tkj laaftmotlon of (k-j). When itaalf la a matrix, than Tv la 
oallad a blook-ToapHta matrix. Itt thia saotloo our attantlon will ba raatrlotad 
to tha oaaa whan tha blook-ToapUta la real, aymmatrlo, peraymmatrlo and poal- 
tlva dafluite. Tha derlvatlona ara vary almllar with the almpla Toq>llta oaaa, 
dw only dlffarenoa being tha raplaoamant of aoalara by matrloaa. Thua wa 
ganarallaa tha algorlttim derived In the prevloua aaotlon. A oonq>lete derivation 
for the Invaralon of poaltlve definite blod(>Toapllta matrloaa oan be found In 
Kutlkov (1966). 

Blodc-Toepllts matrloaa appaar qnlte often In time aarlaa analyala. The 
Blmple Toapllta forma uaually arlaa In one-dlmena tonal Wiener filtering, whera- 
aa blook-Toapllts forma arlaa In two-dlmenalonal filtering. Suppoaa xt , a - 
0,1,. ..,p are random vaotora, each vector having length of (q*fl), then 
the autocovarianoe matrix (atatlonarlty alao aaaumed) la of blook-Toeplltz form 
with dlmena tonality [(p-fl) (q4>l), (p<fl) (q<f 1)] . Weoanaeethat (q+l) la 
the dlmenalon of the aquare aubblooka and ( p**- 1 ) la the number of random 
veotora. 


8.4 Inveralou of Blook-Toepllta Matrtoea 
Development of die Algorithm 

The algorithm given for the almple Toeplltz oaae is generalized here aa 
followB: 

Simple Toeplltz Caae Block Toeplltz Caee 
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Now we can write the iterative solution algorithm for the inversion of the 
block Toeiditz matrices corresponding to the eqns* (3.18) and (3. 19) of simple 
Toeplltz inversione 
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( 1 ) 


( 2 ) 


B,. ■ I 


■ "Otltt.ttP 


Bi “ s**0, 1,». 

..(P-l) 

Aff" §“01 If a a 

• *P 

• • Ei Am § * Of If a a a I (P*1) 


Bi*i,t " Bi,t-x + Bt*i,oBi,i-t X 

a - 1,2, 

Ba^Xfist^x * Bi,p.t ^ Bp«x,o / 

t - 1,2, 

“ App t"0, l,*..,p 



(3.21) 


[ from lymmetry and peraymmetry] 


- iCt 
• b!i 

" IWt 


(3.22) 


i(p-fl)/2 If p odd 

a - (PI+D.....P. where p/g If p even 


Rt-i,i-x ■ “ Ri-x»t-x 


I [ from aymmetry and peraymmetry ] 




Similar to the aimple Toeplitz oaae the aaving in computational time and in 
computer atorage ia very aignifioant. We need only atore one block-row of the 
T* and one block-row of the inverae (Tn* )*’^ tax oaaea where we are not intereated 
in the inverae itaelf but rather in the aolution. A computer program for the aolution 
of (3.23) written in FORTRAN language ia given in Appendix 3.B. 

The computational time by dxia algorithm ia proportional to p*(q-f 1)’ com- 
pared to the {^(q't‘1)^ of the olaeBical method of bordering. Some teat calcula- 
tions on the Amdahl-470 yielded the results :'iven in Table 3. 4 in the solution of 
the linear equation 


T • F - X, N + 1 - (q+l)*(p+l) (3.23) 

where liS is the blook-Toeplitz matrix of dimensionality [(n<fl),(N-i'l)] consisting 
of [ ( q+1 )» ( q-i-l ) ] subblocks, F is the solution vector of length ( N 1 ) and X is 
the dhta vector of length ( N<fl) . Thus for the solution vector we have. 
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(9.84) 


P - (TW*)"‘ X 


IW>1«3.4. CPU Tlmt ud Cor* 8tort|t for tlM bivtrtlon 


of Biook-To^lite Matriooi. 


E 

dimension 

CPU Tims 

(SSC.) 

Storage 

(K) 

P+1 

Q+1 

N+1 

25 

2 

Rl 

MSM 

3 

100 

2 



10 

mllm 





80 




40 

20 



48.88 

93 

Em 

2 

800 

49.68 

40 


We will now ehow that the ooverfanoe matrix of equally epaoed gridded data 
la of block-Toeplita form. If data are gridded and obaervatlons on each profile are 
equally spaced, then distances between profiles do not have to be equal to each 
other In order to obtain a block-Toeplltz form. As an example let us consider the 
three profiles shown in Figure 3. 3. The covariance matrix of this model is of 
block-Toepllts form and is given below. 
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Figure 3. 3. Three parallel Profiles with Equally Spaced 
and Gridded Observations. 
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C(0) C(ii) 

Cfii+H) 

C(0> 

C(fs) 

qrmmetrio 

C(01 

C(a) C(/si* + f*') 
C(i) 

C(^+Xe» ♦ Hf ) 

^nunetrio 

€(•) 
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C(/4s" ♦ 1 

C(2a) 

symmetric 

C(2s) 

rc(3s) Cf/SFreT*) 

C (/9s* + s,*) 

C(3s) 

wymmetrio 

C(3s) 


bi order to elaborate the applloabiltty of the blook-Teopllts Inversion In 
practice, It waa uaed ki the recovery of point and (l^x mean anomalies from 
geold helots. Point anomalies predicted from grldded geo Id heights, which are 
created from the adjusted satellite altimeter data, are shown In Tables 3. 5 and 
3.6. The computation points of Table 3.5 are located on 90-E Ridge, and those 
of Table 3.6 are located In flie Bovin Trench Area. These locations have been 
selected due to liielr large anonuily variations. (l*’x 1*) mean anomalies predicted 
from grldded geold heights are shown In Tables 3.7, 3.8, and 3.9. The (l°x I**) 
computation blocks of the Tables 3.7 and 3.8 are located In the Philippines Trench, 
and those of Table 3. 9 are located In die Calibration area. 

QANOM and CANOM are common In all the Tables and have the following 
meanings I 

GANOM: Gravity anomaly computed from GEM 9 potential coefficients, 
CANOM: Gravity anomaly computed from original adjusted satellite 
altimeter data using the rigorous least-squares collocation. 
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W« ••• from Toblo S. 6 tint tbo RMB dlfbrMotf bttwom tho grtvlty anomalloi 
oomfutid from (ho ortf ImU adjuotod alltmttor data aad tboaa from tlm grtddad data 
(prodlomd) ara 8*0 and 3.0 mgala roapaotivolyt whaa uaing (0 x 85) aad (10 x 40) 
grlda* Thaao RMB dlffaroaooa art wall balow tha ataadard daviattoaa of dm pro- 
dlotloaa thaamaWaa If wa ooaaldar tha ataadard davlatloaa of pradlotioaa to ba 
87 mgala aa bafora. Howavart dm RI0 dlffaraaoaa la Tabla 3.0 ara 27. 6, 18.9, 
aad 17. 8 aigala la dm oaaa (8 x 73), <10 x 38), aad (18 x 20) grlddad data, 
roapaodvaly. Oaly oaa of dm RMS diffaraaoao la alli^dy hlglmr thaa tha 27 mgala 
ataadard davlatloa of pradlotioaa ooaaldarad bafora. 

Wa alaooomputad (l^x 1^) maan anomallaa of blooka 710,711, aad 343 
givaa In Tab’aa 3.7, 3.8, and 3.9 raapaotivaly. Tha RMB dlffaraaoaa batwaan 
nmaa aaornaima oomputad from original altlmatar data and (10 x 40) grlddad 
data ara 8.9, 6.2, 2.0 mgala raapaotivaly. The ataadard davlatloa of (l^x 1^) 
maan anomaly pr^lotlon la about 7 mgala. So tha RMS dlffaranoea above ara 
amallar than the atandard deviation. 

From tha raaulta above, we can oooolude that. In the oaae of altimeter 
data, the RMS dllfarenoaa between pradiotlona from original data and grldded 
data ara generally lean than tha atandard deviation of predlotlona. Thua, we 
oaa oraate a grldded data and predict algnala by tiie fhat Toeplltz algorlthma 
daacrlbad In thla chapter. 
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Ttbto 3.6. Point AaoatllM from OEM 6 Coofflolonto. Origlnol Goold Holghto ond Ortdded 
Ooold Bilifiti. 
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RMS Dlfforoooo Botwoon CANOM and TANOMl la 8.0 mpila 
RMS Dlfforenoa Batwaon CANOM and TANOM2 la 3.0 mfala 


whaia, TANOMl t tha point anomaly oomputad from (6 x 66) grtdded gaoid 
hal^ta. 

TANOM2 1 tha point anomaly oomputad from (10 x 40) griddad geoid 
belghta. 


The grid llmita for pradlotlon deaorlbed In Ihia table arot 
■ -8®, - -14®, Xt ■ 92®, Am " 88^ 
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Table 3. 6. Point Anomalies from GEM 9 Coefficients^ Original Geold Heis^ts and Grldded 
Geo Id Heights. 
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RMS Difference Between CANOM and TANOMl Is 27. 8 mgals 
RMS Difference Between CANOM and TANOM2 Is 16.9 mgals 
RMS Differe:ace Between CANOM and TANOM3 is 17. 8 mgals 


where, TANOMl ; the point anomaly computed from (5 x 72) grldded geold 
heights. 

TANOM2 : the point anomaly computed from (10 x 36) grldded geold 
heights. 


TANOM3 ; the point anomaly computed from (18 x 20) grldded geold 
heights. 


The grid limits for prediction described in this table are; 


<Pn = 30°, <0s = 26°, Xe = 145?75, Xw = 139?75 
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Table 3.7 
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. (l°x 1^) Mean Anomalies from GEM 9 Coefficients, Original Geold Helots and Gridded 
Geoid Heights. 
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R"1S DIFFERENCE BETWEE7 CAROM AN3 T/\nOM IS 5.9 MCALS 

where, TANOM: the (l**x 1®) mean anomaly computed from (10 x 50) gridded 
geoid heists. 

The grid limits for prediction described in this table are: 
iPN = 11®, = 4®, Ac = 128®, A« = 121® 
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Table 3.8. (l®x 1®) Mean Anomalies from GEM 9 Coefficients, Original Geoid Heights and Gridded 
Geold Heights. 
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WS DIFFERE.IC;. BETWEE" C^UtOM /\N TAMOM IS 6.2 HGALS 

where, TANOM : the (l®x 1®) mean anomaly computed from (10 x 40) gridded 
geoid heists. 

The grid limits for prediction described in this table are: 

= 11®, <Ps = 4°, Af = 138®, Ah = 131® 
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Table 3.9 
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R!"S DIFFER-r.CE BETWEE'I CAROM AR9 TAWOM IS 2.6 MCALS 


where, TANOM : the (l®x 1®) mean anomaly predicted from (10 x 40) gridded 
geoid heists. 


The grid limits for prediction described in this table are 
0n=4O?5, (0s=34?5, A£ = 297?5, X« = 291?5 


4. Fourier Tranafonng and Fregueney Domain Collocatiop 


Satellite altimeters such as those on Geos-3 and Seasat-1 satellites intro- 
duced in the last several years have provided multitudes of mearurements of 
gu>id heights. In Chapter 2 we analyzed some of this Geos-3 altimeter data, hi 
the near future other measurements such as gravity gradients by satellite borne 
gradiometers, range and/or range-rate measurements by satellite-to-satellite 
tracking (SST) techniques may be available in large quantities. The processing 
and interpretation of this data by conventional methods is very costly and time 
consuming, even sometimes being impossible. The analysis of such data usually 
involves convolutiou in the space domain, e. g. conventional collocation, but only 
simple multiplications in the frequency domain (see sections 1. 7. 2. 1 and 1. 7. 2. 2 
for convolution) as we will show later in thin chapter. Thus, frequency domain 
methods can be considered as efficient tools in our computations aiad analyses. 
We compute frequency domain representations via Fourier transfoiins of space 
domain representations. 


Fourier transforms can be used for many different purposes. Jordan (1978) 
shows how to use Fourier transforms for upward v'iontinuation, for computing 
anomalies from geoid heights and for solving Stokes' integral, etc. as an appli- 
cation of Wiener filtering (Moritz, 1967), which is a special case of our approach 
to be explained in this chapter. Fourier transforms will be applied on the expres- 
sions in the space domain. Therefore, the theory of least-squares collocation in 
the space domain will be briefly ejq)lained 'oelow just for the purpose of complete- 
ness. For details readers are referred to Moritz (1975). 

4. 1 Least-Squares Collocation 

Least-squares collocation is a method utilizing minimum variance estimation 
for a model of the form (Moritz, 1975, p. 7) 

X = AX + s’ +n (4.1) 


where 


X ; the N-vector of observations (measurements) 
X ; the u-v€ctor of parameters or unknowns 
n : the N-vector of the measuring errors (noise) 
A : the design matrix of dimension ( N x u ) 
s' : the N-vector of signals measured 


The vectors s' and n are purely random and they have zero expectation (average 
or mean value), i.e. 


E{s'} = E{n} = 0 


(4.2) 


Let "s” denote the signal vector of length P to be predicted defined as 
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(4.3) 


,T 

8 » •f) 

and let U8 aaaume zero expectation for signal s and no correlation between the 
signal measured or ttie signal to be predicted and the noise, i.e. 

E{s) * 0, Efs'*n^} * 0, E{e»n^1 « 0 (4.4) 

*ben we can proceed as follows. 

Eqn. (4. 1) can be written in the following way 

X - AXMO (4.5t 

where 0 denotes (NxP) zero matrix and I denotes (NxN) unit matrix. 

It is convenient to substitute 

V « r® 1 (4.6) 

Ls' +n J 

playing the role of "residuals", and 

B » to n (4.7) 

so tttat eqn. (4. 5) can be written us 

X = A X + B V (4. 8) 

Eqn. (4. 8) has the form of condition equaticuswith parameters (Uotila, 1967). 
Here, for least-squares adjustment, we minimize 

V = minimum (4.9) 

where Q is the covariance matrix of the v-vector and may be defined as a 
partitioned matrix 


Q « TCm Cm* *1 (4.10) 

LCg* • Ci*g* C|ja -1 

where Cm = cov (s, s), CsS» = cov (s, s') 

C,», = cov (s', B), Cna “ COV (n,n) 

C,t,*= COV (s', s') 

The covariance between two functions, say f and g , is defined as 

cov(f,g) = E{[f-E(f)l [g-E(g)lM (4.11) 

Thus using the definition above and recalling the aBSunq>tions given by eqns. (4. 2) 
and (4.4), we obtain 
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£(■»•» n. C„ - EfnnM 
Ef«' ■} 


(4.12) 


Cm = E{b8^], C«'*» » 

Cm! “ Efs S*T» Ci*i “ 

The minimization of eqn (4.9) ia equivalent to 

# » i Q"'v - (A X + Bv - X) (4.13) 

where K is the column vector of correlates or Lagrange multipliers. In order to 
minimize eqn (4.13) we proceed as follows: 

|J. = Q”^ - B = 0 

9v 

d4 r 

H A » 0 (4.14) 

AX + ^ -X - 0 

The solution of the equations above : ieid: 

X = [A^ (BQB^)"^ A)"^ A^BQB^"^ x 
v = QBMBQB^"^ (X - AX) 

Substitute eqns. (4. 7) and (4. 10) above for B and Q respectively to obtain 



X = (A^C'^A)"^ A^C"^ X 

(4.15) 


ran rc„» C“^(x-ax)*i 
^ Ls' + nJ L X - A X J 

(4.16) 

where 

5 — Ci>t' + Cm 

(4.17) 


By using the definition (4.11), with eqns. (4.2) and (4.4), it can be easily seen 


C*K — C — C»*,* Cnn 
Cgx “ C*ff , Cxi “ Ci» i 
so we can write for s given by eqn. (4. 16) 

8 = C.x (X- AX) 

The errors of estimations are (Moritz, 1975, pp. 32-33) 

Exx = (A^CxxA)-' 

Egi ~ Cm — Cix Cxx Cxi "f" HA Exx A^H^ 


(4.18) 


(4. 19) 


(4. 20) 
(4.21) 


- 96 - 


’.viiere H » C,» C 

If we auume there are no parameters (i. e. X * 0 ) in the adjustment, as 
we assumed In the second and third chapters in the determination of gravity- 
anomalies and undulations from Geos-3 altimeter data, then we can write (c. f. 
VI. (2.4)) for the signal 


s “ Ci* CxK a (4. 22) 

and for the errors on signal 

Eii * C(f “ Ctx C*a Cm (4, 23) 


The equations given in this section solve our problem in the space domain. 
The inversion of the covariance matrix of observations , which has a dimen- 
sion equal to the number of observations ( N) , is very costly and time consuming 
for large N . Therefore, we seek another solution, which is faster and c eapcr 
than the space domain solution, i.e. the frequency -domain solution. It io this 
alternative which is to be discussed in the next sections. 


4. 2 Frequency Domain Collocation 

hi the previous section we derived the space domain equations of least- 
squares collocation, b what follows, we will derive the corresponding equations 
in the frequency domain. Thus the least-squares collocation is carried out in the 
frequency domain and solved for the desired quantities. The application of the 
inverse transform (Fourier) yields the corresponding desired quantities in the 
space domain. All the required transforms are carried out by the fast Fourier 
transform (FFT) to facilitate and speed up the computations. 

We have defined the discrete Fourier transform (DFT) pairs by eqns. (1. 65) 
and (1.66) in section 1. 5. Here we use a slightly different definition for the DFT 
pairs ejq)ressed as follows; 


N-l 


X„ = -^ x(k) exp(-i2rrkn/N), (direct Fourier transform) (4.24) 

n = 0,1,...,(N-1) 


. N-l 

x(k) = Yt Xa exp( i2TTkn/N), (inverse Fourier transform) (4.25) 

k =0,1,. ...(N-l) 

The equations above are derived from eqns. (1.65) and (1.66), first by multiplying 
both sides of (1. 65) by /IT in order to obtain 


Yn = Xn' = ^^^J^^x(k) e;q>(-i2TTkn/N) 


(4.26) 
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and then subatitutlng X.' * Ya//IT in eqn. (1.66) to obtain 

, N-X 

Eqns. (4. 26) and (4. 27) are equivalent to (4. 24) and (4. 25) respectively, and can be 
soen following the replacement of Y. by X, above. 

Eqn. (4.24) can be written in a matrix form as 


X » Aj'X 

where x » ( Xo, Xi,...,X<-ir, x = (x(0), x(l), . . . , x(N-l) 

g-snjoo/« ^-aTTixo/i g-an i(n-i)o/n 

„-anjoi/N -anJuA „-9 tti(n-i)x/n 

: : : 

gSnio(«-i)/« -aTrix(N-i)/N g-ani(N-i)(N-i)/N 


1 

7W 


(4. 28) 


(4. 29) 


Here Ai is the transform matrix and equal to the Fourier matrix defined by 

F = [ Ffcnl = -;— [ exp(-2n lkn/N)l, 0 ^ k,n ^ (N-1) (4. 30) 

/W 

io we can write Ai = F (4* 31) 

As we mentioned before, we used the space domain equations (4. 22) and (4. 23) for 
the determination of gravity anomalies and their errors from geoid heights. Now 
we can derive the corresponding equations in the frequency domain. The Fourier 
transform pairs for observations, signals and covariance matrices will be denoted 
as follows. 


X < — > X, 

Cix < '■”> Csx* 


3 <— >> S, 
Cxx < — > Cxx I 


Cxf < ^ Cx8 

Eat ^ ^ E33 


(4.32) 


The vector of measurements (x) and the vector if signals ( s) are trans- 
formed into the frequency domain throu^ 

X ~ A^ X , S “ A]_ s (4. 33) 

where the transform matrix Ai is defined by eqn. (4.29) being equal to the 
Fourier matrix defined by (4. 30). In order to obtain the frequency representations 
of the covariance and cross-covariance matrices Cxx» C,* $ Cix, we recall the 
definition of covariances and cross-covariances given by eqn. (4.11) 

Cxx = E{[X-E(X)1 [X-E(X)1+} (4.34) 
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where (t) denotes the complex conjugate and transpose operation, 
eqn. (4. 28) for X in (4. 34) to obtain 

Substitute 

Cxx - E([AiX - E(AjX)l (AxX - E(AxX)l'^] 


» Ai Ef(x - E(x)l [X - E(x)J } Ax"^ 


Since the expectation of the terms inside the braces is equal to C 
covariance matrix cl observations, we can write 

XX • the 

Cxx “ Ax Cxx Ax 

(4. 35) 

Similarly, we can prove that 


Csx = Ax C*x Ax*^ 

(4. 36) 

C)(s ~ Ax Cxi Ax*^ 

(4.37) 

Css “ Ax Cfi Ax*^ 

(4.38) 

Ess “ Ax Eli Ax"^ 

(4, 39) 


Having transformed the space domain quantities into the frequency domain, 
we can write the frequency domain collocation solutions corresponding to eqns. 
(4.22) and (4.23) of the space domain solutions 

S = Csx Cxx X (4.40) 

Ess “ Css - Csx Cxx Cj(s (4.41) 


Following the computation of the signal vector S and its error Ess in the 
frequency domain, we can find the signal vector s and its error Em in the space 
domain by applying the inverse transform as explained below. From eqn. (4. 33), 
we can write 

s = Ax^S 

Since Ax = F is an orthogonal matrix (recall the orthogonality relations of ex- 
ponential hinctions as explained in the first chapter), we have for the signal 

s = a| S (4.42) 

and similarly for the error matrix 

E„ = aJ* Ess Ax (4.43) 
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1 


1 


1 


1 7^ T T 




T 




We see from eqos. (4.40) snd (4.41) that we still have to invert the matrix 
CxK » which has dimension (N x N)» N being equal to the total number of obser- 
vations. We transformed everything into the frequency domain to avoid the 
inversion of the co/ariance matrix Cx» $ which also luu dimension (N x N). The 
computations of Cxx* Cix* and Cx« from eqns. (4.35), (4. 36) and (4. 37) respec- 
tively, require the product of three (N x N) ftill matrices in each case. So if we 
go blindly, we will have to perform mai^ more computations compared to the 
direct collocation solution in the space domain. However, we will show later 
that a frequency domain matrix (FDM) is diagonal If the corresponding space 
domain matrbc (SDM) is circular, and asymptotically diagonal if the SDM is 
of the Toeplitz form. In addition, we will show tha^ uiagonal c.’ements of 
the FDM is the Fourier transform of the first column of corresponding SDM. 
Hence, we do not have to compute the FDM by the dii4>« t transfoi.uiatioo of the 
corresponding SDM, for example, as ineqn. (4.35), which requires the product 
of three (N x N) matrices. 


4. 3 Diagonalization of Toeplitz and Circulant Matrices 


Suppose we have a stationary covariance function and equally spaced (or 
sampled) observations along a profile, then, as we explained in the third chapter, 
the covariance matrix of observations are of the Toeplitz form. In addition, if 
this profile forms a complete circle, then the resulting covariance matrix is 
circular. /•. circulant matrix, say Tc , is one having the form 


Tc 


to k 

tw-l to 
tw-a tn-i 

• m 


ta 

ti 

to 



(4.44) 


ti 


»a 


t. 


to 


The circulant matrix Tc is a special type of the simple Toeplitz matrix Tn_i 
given by eqn. (3. 1) such that t_n =• tN_k , k = 1, . . . ,(N-1). The diagonalization 
of Tc and Tn_i will be given in order to see that the eigenvalues of Tc and Tn_i 
( asymptotically) are the diagonal elements of their frequency domain repiesentaticns. 


4. 3. 1 Diagonalization of loeplitz MStrices 

As in the third chapter, we denote a simple Toeplitz matrix by Tn-i . The 
dimension of Tn_i Is (N x N) and defined as in eqn. (3. 1) such that its elements 
on symmetric diagonals are identical and its diagonal elements are equal. Asrume 
the covariance function t^ of Tn_x is absolutely summable, i.e., 

00 

E kl ^ ” (4.45) 

k=-oo 

then, for large N, Fuller (1976, pp. 133-138) proves that 
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(4.46) 


— TT^ 
tty 




1 F ■ Diftg (Ao» Aj,..., Am-i) 

whore F is defined by eqn. (4.30) and V. are the eigenvalues of Tn. The 
eiger values (A#) are defined as follow;*; 

Ab - 2tt t(o)»), oJb « 2rrn/N, n ' 0, 1, . . . , (N -1) 

where t(cOg) is the spectral density of the t^ e^qpressed as 

■■= fr Z exp(-iw„k), n » 0,1,...,(N-1) (4.47) 

* tt kat.OB 

Thus, we can write for eigenvalues 

00 

An ® n = 0,l,...,(N“l) (4. 48) 


Since we cannot carry the summation in eqn. (4.48) to infinity, and becasue 
t^ is absolutely summable, i.e. t|, gets smaller with increasing |k|, we will 
assume: 


(1) ti, = 0 for N^: we can write for the eigenvalues of Tn_i 

= to + X (1, 

Hence X = /IT^F^'t* 


(4.49) 


where A [Ao,Ai , . . . ,Xn_i ] * ^ ~ [^»(^ 3 )* • • • *(^-i'^^i)]* and 

F is the Fourier matrix defined by eqn. (4.30). 

In the following sections, we will come across products such as FTn-i f"^ 
instead of F'^Tn_iF as given in eqn. (4.46). Therefore, we will define the 
vector of the eigenvalues of FTn_i F"^ as follows: 


A = /FrF*^t‘ 


(4.50) 


(2) tu = 0 for k > m, where m ^ N/2 : with this assumption, Tn_i becomes 
(to be denoted as Tm_j ) 


Tn-1 


t 

to 




t_, 0 

t> «....t^ to t_^ .... t^B 
• • . 

o'-.. 


1. 1.,1 


‘t- 

I 


(4.51) 
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If we fill the zero elements of Tn-j In such a way that the resulting matrix, 
say Tc , is clrculant, 1. e. 


where 



t - 

0 . 

t h ^ 


Co 

C. .... Cn_i 

•o X • • • • 

• 

t i'— 1 

t-.*’ 

* Q • • • • tg 


9n-i 

C 0 • • • • C^ — 3 

« 

t. 

to 

t_l . 

• 




0 ., 

: • • . 


to . 

...... .t_„ 




: 

t 

•0 

t . .. 

to ' 


Cl 

C« • • • • Cq 

*'-1 • • • 







/ t_K for k = 0,1,. . . ,m 
Ck = \ tfj-k for k = (N-m), . . . ,(N-1) 
( 0 otherwise 


(4. 52) 


(4. 53) 


then Tn_i and Tc are asymptotically equivalent which can be proven aa follows 
(Gray, li>77, pp. 27-29); 

■ 

Urn |Tn_i-TcP = lim^ k [ + \u\- ] 

^ lim ^ t [ |t,|^ 4 - |t_J®] = 0 

from the sumniabil.ty of covariance functions. Thus 

lim (Tn^ -Tcl - 0, i.e lim Tn - Tc (4.54) 

N -♦cT' N “ 00 


In what follows, we will see that the diagonal iKation of a circilar matrix is 
simple and its eigenvalues can be found easily. Since is asymptotically 

equivalent to the circular matrix Tc , the eigenvalues of Tn_x are equal to those 
of Tc and Tn_x can be diagonalized similar to Tc. 


4. 3. 2 Diagonalization of Circular Matrices 

We defined a clrculant matrix Tc by eqn. (4.44). The eigenvalues, X„, 
and eigenvectors, Yq, of Tc satisfy the following equation. 


Tc * Yn = • Y„ , 

n = 0,1,...,(N-1) (4.55) 

where Y„ = [yin. Ysn. • • • .y(N-i) n 1 
that 

_ _ p-lSnn/N 

Let r„ be a root of the equation = 1 such 

n = 0,1,...,(N-1) (4. 56. a) 

and let us set 


Vin = r„^ 

j,n = 0,J :N-1) (LCr-.b) 
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then eqn. (4. 55) can be written 


Tc 


r° 

r„^ 


.N-l 


= x„ 


r„‘ 


.H-n 


(4. 57) 


Now substitute eqn. (4. 44) for Tc to obtain 


Co + Cl r„^ + 
CN-ir„° + r*,r„^ + 


+ CN-irn 


N-i _ 


= X„rn“ 


N-X X 1 

+ Cm 1*11 "■ An Tn 


(4. 58) 


Cl r„° + Ca r„^ + + co rU-^ = X„rn 


-N-* 1 


and multiply the first equation of (4. 58) by Tn , the second by and so forth. 

Then use = rn in order to see the equality for each equation if 


N-l 


x„ = E 

k«0 


Ck • To 


The substitution of eqn. (4. 56. a) in (4. 59) yields; 

N-l 

Xn = Ck e-‘»’^kn/N, n = 0,1...., (N-l) 


(4.59) 


(4.60) 


Hence, the characteristic vectors are given as 

Yb = [i,3-i2nn/N^ e-isn(N-i)n/N ]T^ n = 0,1,..., (N-l) (4.61) 


and the orlhonormal matrix, which also diagonalizes the original matrix Tc , of 
eigenvectors is f,. ->n as 


F 


= ^[Yo. Y. 


Yn_i1 


(4. 62) 


Eqn. (4.62) is equivalent to eqn. (4. 30) defined previously, i.e. the Fourier matrix. 
Since F is the orthonormal matrix of eigenvectors, we can write (Byron et. al. , 
1969, pp. 120-124) 


f"^ Tc F — Biag ( Ao » Xi , . . . , Xj ^1 ) 

Hence, 

^ — [ Xo , Xi , . . . , — 1 J — /N' * F * c 


(4. 63. a) 
(4. 63. b) 


•I* 

Here F denotes the complex c onjugate transpose of F as usual and c is the first 
row (or column) of Tc. Similarly, the diagonal elements of FTc F^ are expressed 
as follows 


X = /F F+ c 


(4. 63. c) 
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Thus wo have proved that the Fourier matrix ( F) given by eqn. (4.30) Is 
the orthonormal diagonalizing matrix formed by the eigenvectors of elt' nr a 
Toeplltz matrix or a circular matrix. The eigenvalues can be computed simply 
by eqn. (4.50) for Toeplltz matrices and by eqn. (4.60) for circular matrices. 


4.4 Frequency Domain Collocation on 0>ie-DImenslonal Region 
4.4.1 Large N and the Classical Wiener Filtering 


In section 4.2 we derived the least-squares collocation expressions in the 
frequency domain, and In the previous section we showed that the transform 
matrix Aj = F defined by eqn. (4.29) is, in fact, the diagonalizing matrix of 
the covariance matrix Cxx» Cx» > C,x. Therefore, we can write for the frequency 
domain collocation covariance matrices Cxx» Cxs. Csx. and Css . 

XX XX xx 

Cxx “ Diag ( A.0 f \\ 9 • • • I ^—1 ) 


Xt x» x« 

Cx3 ~ Dia.g (Xof »•••• A^~i) 

Csx “ Diag (Xo , X 1 , . . • , Xw_i ) 


(4.64) 


Css ” Diag ( Xo » Xi , . . . , 



where X”, X”, X**, X“, n = 0, 1, . . . ,(N-1) are the eigenvalues of Cxx,Cx«.C«xi 
and C,» respectively. Thus, we can write for the frequency domain collocation 
solutions (expressed by eqns. (4.40) and (4.41)) 

s„ = xr* x„/xr (4.65) 

n = 0,1,...,(N-1) 

(Ess)„ = X*„* - X^X^Vxr (4.66) 


where S„ is the n-th signal of S vector and (Ess)„ is the corresponding error, and 
Xn is the n-th value of X vector, i.e. S - [So,Si , . . . ,Sn_i] \ X = [Xo.Xi , . . . , Xn.J^ 


The eigenvalues in eqn. (4.64) are computed from eqn. (4.50) when tlie cor- 
responding space domain matrix is of the Toeplltz form and from eqn. (4. 63. c) 
when it is a circular matrix. The equations given by (4. 65) and (4. 66) lead us to 
the very well known "the classical Wiener Kolmogorov filter". 


Having found S in the frequency domain we can easily obtain the corresponding 
s-signal vector in the space domain applying the inverse transform as expressed by 
eqn. (4.42). The computation of the error matrix Ess from Ess tli rough eqn. (4.43) 
is also very simple: Consider the diagonalization of the circular matrix Tc by 
eqn. (4. 63. a). Now suppose the diagonal elements of F^Tc F, i.e. the eigenvalues 
X„v are given and the circular matrix Tc is asked. Then the first row c of Tc 
is found from (4.63.b) by multiplying both sides of the equation from left side by 
F^ to obtain 
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c = 1 ^ 

- 7f ~ 

since Esj = FEjs has a similar form with (FTc F"^) in addition to its 
diagonality, the error covariance matrix of signals (Eas) is circulant and its 
first row S, is computed as follows 

S.=-^Fa” (4.67) 

where d” = [do » d^ 1 are the diagonal elements of Eas. The com- 

putation of the remaining elements of E„ is straightforward by recalling the 
definition of circular matrices given by eqn. (4.44). 

Thomas et. al. (1976, pp. 3-4) relates the overall estbnation error to a 
single quantity, namely to the rms estimation error, (7, defined by 

as® - ^ Efs s'] - f E{(Sr-s)''] (4.68) 

N N k»i 

where N is the length of the s-signal vector to be predicted 
St is the vector of true signals. 

Eqn. (4. 68) is equivalent to 

as* - -J- Trace (Ess) (4. 69. a) 

N 

br by the Parsenval's theorem mentioned in section 1.3.4, it is equivalent to 
. '^1 ss 

as* = Trace (Ess) =. Xk (4.69.b) 

N k = 0 

Thus we have completed the solution of the s-signal vector and its error. The 
total number of computations for the frequency domain least-squares collocation 
explained above is proportional to N log ^ N multiplications and additions compared 
to N® in the case of the space domain least-squares collocation. 

Moritz (1967) introduced the frequency domain method explained above as 
"the least-squares filtering", and Sjoberg et al. (1977) applied it for the picdiction 
of mean free-air gravity anomalies from altimeter data along one arc at a time. 

In his computations, Sjoberg assumed the covariance or cross -covariance t,, to be 
equal to zero for k > M (M = 70) compared to t^ = 0 for k N in our solutions. 

Following the derivation of the frequency-domain algorithm, it was necessary 
to test it. In order to do so, we selected two arcs, which are the 11-th and 12-th 
arcs of Figure 2.2. Our aim was to determine free-air point gravity anomalies 
from altimeter data at data points along one arc at a time. We followed the com- 
putational steps as outlined below: 
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(1) The data vector of geo id heights, say x(k), k » 0, 1, . . . ,N-1, was 
transformed Into the frequency domain by eqn. (4. 33). 

(2) The covariance matrix Cxx of geoid heights and the cross -covariance 
matrix C,x between observations (geoid heights) which were computed by using 
the empirical covariance functions given by the subroutine COVA of Tschemlng 
and Rapp (1974) and signals were transformed into the frequency domain by 
eqn. (4.50). 

(3) The frequency domain representation of the signal vector was computed 
by eqn, (4.65). 

(4) Finally, the desired s-slgnai . in the space-domain was computed 
by eqn. (4.42). 


For our examinations, we first considered a 682 km long segment of the 
11-th arc and denoted it ARC 11-A, and a 1380 km long segment of the 12-th arc 
and denoted it ARC 12 -A. Secondly, we considered Arcs 11 and 12 in their total 
lengths. These arcs and the statistics between the space domain collocation 
solution (SDCS) and the frequency domain collocation solution (FDCS) are described 
in Table 4.1. The free-air point anomalies determined by the SDCS and FDCS 
are shown in Figures 4.1, 4.2, 4.3, and 4.4 for the arcs 11-A, 12-A, 11, and 
12 resp>ectively. Table 4. 1 and Figures 4. 1 to 4.4 do not include about 1. 5% of 
predictions made at the beginning and at the end of each arc. We deleted them due 
to their large differences from those of SDCS. This big difference Is caused by 
the negligence of off-diagonal terms in FDCS known as "edge effects" (Thomas et. 
al. , 1976). 


Table 4. 1. The Statistics Between SDCS and FDCS Along 
The Arcs 11-A, 12-A, 11, and 12. 


ARC 

No. 

Latitude 
<Pi. (Deg.) 

Longitude 
Xi , Ag (Deg.) 

Length 

(km) 

No. of 
Obsns. 

RMS 

Diff. 

(mgals) 

Mean 

Diff. 

(FDCS- 

SDCS) 

Max. 

Diff. 

(FDCS- 

SDCS) 

11-A 

12.848, 18.197 

309.957, 306.844 

681.9 

50 

1.8 

n 

- 4.3 

12-A 

-0.024, 10.881 

146.530, 140.565 

1380.2 

100 

11.7 


-42.3 

11 

12.848, 39.751 

309.957, 291.194 

3511.2 

254 

2.6 


8.4 

12 

-0.024, 32.466 

14S.530, 126.823 

■ . , 1 

4162.3 

300 

5.3 

-0.6 

-33.6 
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Figure 4. 3. Free-^lr point anomalies from SDCS and FDCS along the Arc 11. 
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Figure 4.4, Free-alr point anomalies from SDCS and FDCS along the Arc 12. 



4. 4. 2 Small N and the Windowed Frequency Domain Solution 


In the preceeding section, we asc ^ned a long data length (large N) over 
an interval longer than correlation distances and used asymptotically equivalent 
expressions for the eigenvalues of the space domain covariance and cross-co- 
variance matrices in order to compute the corresponding asymptotically diagonal 
frequency domain covariance and cross -covariance matrices, hi the case of a 
short data length and small N the asymptotic equalities mentioned above are 
no longer valid. Therefore, in this section, we will introduce the windowed 
frequency domain collocation algorithm. This algorithm retains the necessary 
computational speed to solve large problems but uses data given over an interval 
shorter than correlation distances. The frequency domain representations of 
covariance and cross -covariance matrices become band-diagonal here and 
diagonal in the limit as N gets larger and larger. By "windows" we try to 
minimize the number of super-diagonal bands, which is proportional to the side 
lobe energy (for the frequency domain representations of windows see section 
1.7.3. 1) of the data spectrum (Heller et al. , 1977, p. 13). In what follows, we 
transform the "windowed data" into the frequency domain and perform our 
solutions there in order to compute the s-signal vector and the En-error 
covariance matrix. The solution of the s-signal vector given by eqn. (4.22) is 
achieved in two steps ; 

- the solution of y = x 

- the solution of s = C,x x = C«x y 


4.4.2. 1. The Computation of y = Cxx x 

We want to solve a system of equations 

Cxx • y = X (4. 70) 

for the y-vector. Since Cxx is usually a full matrix, this solution would require 
approximately N® computer operations. In order to reduce the number of com- 
puter operations, we first multiply the data vector, x, by an appropriate window 
matrix to control the sidelobe energy of the spectrum of the data vector, then 
transform it into the frequency domain, i.e. 

X = Aa • X (4.71.a) 

and for the frequency domain representation, Y , of the solution vector, y 

Y = (Aa+)"^ -y (4.71.b) 

where Aa is the transform matrix defined as 

Ag = F w (4.72) 
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Here F is defined in eqn. (4.30) and w is a window matrix (usually diagonal). 
For the covariance matrix of X (analogous to the derivation of eqn. (4. 35)) we 
can write 


'XX 


— Aq CxX A< 


= F w Cxx w F^ 

Hence eqn. (4. 70) takes the following form in the frequency domain 

Cxx Y = X 


(4.73) 


(4.74) 


We solve the system of linear equation above for Y , and we compute the y-vector 
in the space domain from eqn. (4.71.b), 


y = aJ Y 


“ w 


Ft Y 


(4. 75) 


Why solve the system of linear equations in the frequency domain such as 
defined by eqn. (4. 74) instead of solving the corresponding equations in the space 
domain such as defined by eqn. (4.70)? The reason is that Cxx is nearly band- 
diagonal and the solution of eqn. (4. 74) can be carried out inexpensively by an 
efficient algorithm such as the banded-C hole sky decomposition (Forsythe et ai. , 
1967). 

In order to show the almost band-diagonal form of Cxx we write eqn. (4.73) 
in the following form (Heller et al. , 1977, pp. 55-59) 


Cxx “ 2 K Fsn w Cxx w Fsn K 


(4. 76) 


where: K is the sampling matrix defined as 


K = 


10 0 . 

0 0 1 . 


.0 0 0 

.0 0 0 


(4. 77. a) 


L? ® ® 0 1 Oj (N X 2N) 

Fa N is (2N x 2N) Fourier matrix defined as 

FaN = [(F3N)ta] "7^ 0 k,n ^ (2N-1) (4.77.b) 

w is the (2N x 2N) extended version of the w window matrix and defined as 

; = ^ (4.77.C) 

lij ^ (2N X 2N) 

Here, w is a diagonal matrix with diagonal elements (wo, Wj , . . . , Wn_i ) 
defined by the window function, and 0 is (N x N) zero matrix. 


-no- 


(4.77.d) 


Cxx is a (2N X 2N) circular matrix defined as 



Here Cix is a Toeplitz matrix with elements 


! Cxxl(je-N)*Atl for jt> 0 
0 for Jt = 0 

Cxx((je+N)*AtJ for Jt<0 


By substituting the definitions of K, FgN t w» and Cxx above in eqn. (4.76), 
we can see that eqn. (4.73) and (4.76) are equivalent. Eqn. (4.76) can also be 
e;qpressed as 

Cxx ~ 2K Fgn w Fsn Cxx Fai Fsn w Faw 


where 


A A A » 

= 2 K WCxx W 

A A A A 

W = FaN W Faw * "^’xx — Faw Cxx Fan 


(4.78) 
(4. 79) 


We have shown in section 4.3.2 that the frequency domain representation of 
a circular matrix is diagonal as given, for example, by eqns. (4. 63. a) and (4.63.b) 
From eqn. (4. 63. a) we can write for the circular matrix Tc 

Tc F • Diag ( Xo»Xx , . . . ,A(g_i ) F^ (4.80) 

A 

If we let N go to 2N,^ then Tc of eqn. (4. 80) and W of eqn. (4. 79) have similar 
forms. Therefore, W is also circulant with the elements 

[ W]itn “ n » 0 ^ k, n ^ (2N~1) (4. 81. a) 

A 

we can write for the first row, say O , of W 

Q “ ( Oo» » • • • » OsN-i)^ ~ (4.81.b) 

as can be derived from eqn. (4. 63.b), where (^) is the vector of the diagonal 
elements of A' ^iven by eqn. (4. 77.c). It also — follows directly from eqns. 

(4. 63.a-c) that Cxx given by eqn. (4. 79) (recall that ^xx is circulant) is 
diagonal as expressed below 

Cxx “ Diag (Wo» » • . « f <OaN-i ) (4. 82. a) 

Actually, the U)„'s are the eigenvalues of Cxx given in eqn. (4. 77. d) and they can be 
computed using eqn. (4. 63. c) 

to = Fat c (4.82.b) 
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where w 
c 


= ( otef coi . . . . , cOgN-i)^ is the power spectrum of 
Is the first column of the circular matrix 6xx • 



and 


A 

Cxx 


Now substitute eqns. (4, 77. a), (4. 81. a), and (4. 82. a) for K, W, and 
respectively, In eqn. (4.78) and use some simplifications to obtain 


aN-i 

(<5xx)jk = ^ j^Q 


0 s- j,k 5^(N-1) 
= OsN-X 


(4.83) 


w^here Ojj and cOjj are defined by eqns . (4. 81. b) and (4. 82. b), respectively, and 
Ql represents the complex conjugate of Ojj. The equation (4.78) Is now a very 
efficient formula for the calculation of Cxx • This comes about because when we 
select a window which has a mainlobe and a series of very small sidelobes, then 
Cxx has an almost band-diagonal structure with some elements In the upper right 
and lower left comers. Hence, the solution of eqn. (4.74) can be carried out by 
the band-diagonal and comer Implementation of the Cholesky decomposition, which 
requires numerical operations proportional to (Nb* • N), where Ne is the width 
of the band, compared to N^ for conventional solutions. 


The out-of-band elements of Cxx can be made very small by selecting an 
appropriate window and matrix bandwidth (usually less than 10) so that we can 
neglect the out-of-band elements In our solutions. If Cxx with zero out-of-band 
elements Is denoted by Cg , then we can state that Cq Is an approximation to 
the exactly transformed matrix Cxx • This approximation affects the solution 
vector y of eqn. (4.70). The magnitude of the error on y can be controlled by 
the choice of window and matrix bandwidth as mentioned above. The approximation 
above is the only error introduced In the calculation of the y-vector. 


There are many windows one can use for tlie purpose of minimizing the apj^ rox- 
Imation affects. Here we will only examine the Kaiser window which has proved to 
be very useful in signal processing applications. The Kaiser window is expressed 
in the space domain as follows: 


where 

and 


Nh 


Wx(n) = Io(^A-(n/N„)'^)/lc(/9), 
f(N-l)/2 if N odd r / ^ _ 

I N/2 if N even’ 


-Nh ^ n ^ Nh 


k = o 



(4. 84) 


j8 is a constant that specifics a frequency response trade-off between the 
peak height of the sidelobe ripples and the width or the energy of the 
mainlobe. 


The space domain and frequency domain representations of the Kaiser window for 
= 2,4,6, and 8 are shown In Figure 4.5. 

Having examined the Kaiser window, we return to the elements of Cxx defined 
by eqn. (4.83). The right-hand side of this equation is the convolution of 0^ with 
the wel^ting function 
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KRISER HtNOOH 








(4. 85) 







Jt-3k 


This weighting function lb very small and consequently negligible when |k - j j > Ni 
and |k • j I < (N - Nb ) which correspond to the product of a malnlobe and sldelobe or 
two sldelobes. It Is for this reason that we have tried to find a window function 
with very small sldelobes. Asa result, we only compute the band-diagonal and 
upper right and lower left comer elements of Cxx shown In Figure 4. 6. Hence, 
we have for the super-diagonal elements 


fCxxlk,k+. = 2 f^A-ak^X-80,+.) 



k (N-m-1) 
m < Nb 


(4. 8C.a) 


and for the upper right corner elements 

a-i 

[Cxxlk.j = 2 f^je-ak^i-3j * 


r 0 < k < Ne 
l(N-NB+k)£ j < N 


(4.86.b) 


The lower diagonal elements of C;;* are just the complex conjugates of the corres- 
ponding symmetrical elements jf the upper diagonal elements. 


If we let w (window matrix) go to I (Identity matrix) and Nb go to zero, 1. e. 
only the elemei ‘s of the main diagonal are retained, then the algorithm presen';ed 
above reduces L "the Generalized Wiener Filtering" described In section 4.4. 1. 



Figure 4.6. The structure of Cxx under Fourier transformation. 
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4. 4. 2 , 2 The Computation of b ■ * V 


For the computation of (C.x* y) we will utilize the fast Toeplitz matrix 
multiplication. Thus, the y-vector will be multiplied by the Toeplitz matrix 
(\x such that this algorithm will take many less operations th;m the standard way 
of multiplying the y-vector by Ci*. Suppose the cross-covariance matrix 
has dimension (M x N), then the matrix-vector multiplication rt>quires roughly 
(M • N) multiplications and additions. If M is very small comparted to N then 
any classical matrix-vector multiplication can be used. However, if M is large 
then the fast algorithm presented below would be more efficient. 


In the preceedlng sections we assumed M ^ N so that Cix is a Toeplitz 
matrix. If M N , then C,x is a partition of a Tot;plltz matrix. In the case of 
M < N we extend C,x in such a way that the resulting matrix, to be denoted as 
Cix, is a square Toeplitz matrix. The first M elements of the product (C»x* y) 
contains the desired product (Cxx* y). 

In order to achieve a fast multiplication of (Cix* y) , we Imbed Cax in a 
(2N X 2N) circulant matrix, say Tc , as explained below. The cross -covariance 
function Vix(k *A0) between signals and measurements is extended to the range 
-(2N-1) ^ kr. (2N-1) by the definition 


E.x(k*A|^l) 


y.x(k*Ai/)) for 0 ^ k (N-1) 

0 for k = N 

%,[(k-2N)*A0] lor (N+1) <k ^ (2N-1) 


(4.87) 


and E«x(-k*A({l) = Eix [(2N-k)‘Ai()] for 1 k ^ (2N-1) 

where A0 is the spacing between s and x. The matrix Tc with the elements 

(Tc)jic = E.x [(k-j)*A0] (4.88) 

is circulant with the property. 

^'(o) ' (“"d'') 


For the fast matrix-vector multiplication, we transform the vector ( q) and 
Tc into the frequency domain as explained in section 4.2 (analogous to eqn. (4. 33) 
and eqn. (4.63) for(^ ) and Ic respectively) to obtain 

Y = [Yo.Yi Y3 n-i]' = (^-90) 

and 

Tc =5 Fqn Tc F^n ~ DUig( Xot Xi t • . . ) 

, ^ . (4.91) 

^ [ ^0 >XT^f*.«»XsN-i] v2 N" Fa N c_ 
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... JUft . 


where Is defined in eqn, (4. 77. b) and C Is the first column of Tc . Hence, 
wo can write _ 

Tc'*Y = = I YoXo, YjXi, . . . , YaN-iXaN-il (4.9i^) 

Finally, we apply the Inverse Fourier transform to obtain 


f-a 0 = 


F^n 


Xo Yo 

Xi Yi 
XaN-iYaN-i 


(4.93) 


The first M elements of the product (C»x* y) is the des Ired product (C«x* y)» 
so multiplication Is completed. The number of operations of the entire procedure 
amounts ot three FFT applications, therefore. It Is proportional to (N loggN). 


4. 4. 2.3 An Application of the Windowed Frequency Domain Collocation 

In order to demonstrate the efficiency of the windowed frequency domain 
collocation algorithm described above, we selected an arc of altimetry (12 th 
arc of Figure 2.2), which is 4162 km long and has 300 observations (N = 300), 
and computed the (5°x 5°) free-air mean gravity anomaly of block 711 (<Pn = 10°, 

<f>i =5®, Xe = 137°, Xw = 132°) shown in Figure 4. 5. We also computed the (1° x 1°) 
anomalies inside 711 from the gcold heights along Arc 12 by the following techniques. 
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\ 

* 


if 


(1) Convention^ L^^-Squajres CoUqc^lon 

Here mean anomalies were computed from eqn. (4. 22) as explained in 
section 4. 1. 


(2) Clas£i^nVi^^ FiM^^ (i.e. Frequency Domain Collocation with I^arge N) 


The y-vector and Cx* matrix of eqn. (4. 22) were transformed Into the 
frequency domain by eqns. (4.33) and (4.50) respectively, to obtain 


Y = C7}X = (Diag(Xo,Xi X 

”lAo 0 0 

0 lAi 0 

0 0 ........ lAti— 1 


(4. 94) 


then an inverse fast Fourier transform was aiq>lied to Y to yield 
y = = f'^'y 


(4. 95) 


Finally, we computed the s-signal vector by 
s = C,x y 

where s represents the signal vector of mean anomalies and Cix represents the 
cross -covariance between mean anomalies and observations (geoid heights). 

(3) Wto^we^Fre5p^]My_I^m^nJ3cl^^ 

The x-vector and Cxx matrix of eqn. (4. 22) was transformed into the frequency 
domain by eqns. (4. 71. a) and (4.73) respectively, to form the following system of 
linear equations 


Cxx Y = X (4.96) 

We solved the linear equaticms above for Y using the Cholesky’s decomposition as |i 

ejq)lained in section 4. 4. 2, then we computed y = C,J^x vector by an inverse 
transform defined by eqn. (4. 75) , i. e. 

y = w F^ Y 

Finally, we computed the s-signal vector of mean anomalies by 

s = Csx y 
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The RMS difference of 25 (l°x 1®) mean anomalies between least-squares 
collocation and Wiener filtering is 0.87 mgals compared to 0. 16 mgals between 
least-squares collocation and windowed frequency domain collocation. However, 
both RMS differences can be considered small con pared to die 15 mgals standard 
deviation of the (l®x 1®) anomaly predictions. Here we have recovered about 3. 5 
mgal anomaly information with respect to the GEM 9 surface. The small magni- 
tude of the recovered anomaly masks the efficiency of the algorithm. To achieve 
a considerably higher magnitude we considered block 614 (<fin - 15®, = 10®, 

At = 139®, Ax = 134®), through which the 12th arc passes. Unfortunately we 
recovered only 5. 2 mgals with respect to the GEM 9 surface, so it did not help 
us much to check the efficiency of the algorithm. 


4.5 Windowed Frequency Domain Collocation with Sine-cosine Coefficients 

In section 4. 4. 2 we have seen lliat the frequency domain covariance matrix 
Cxx of windowed observations has a bnad-diagonal shape with some elements in 
the upper-right and lower-left comers. In addition, the elements of Cxx are 
complex valued, hi order to use the standard banded Cholesky alg orithm for the 
solution of the system of linear equations given by eqn. (4.76), we have to elimi- 
nate the upper-right and lower-left comer elements of Cxx hi addition to real- 
valued band-diagonal elements. If we use a sine-cosine transform (Fuller, 1976, 
pp. 135-137) instead of a Fourier transform, then the resulting new covariance 
matrix Cxx is real-valued and band-diagonal (Heller et. al., 1977, pp. 38-39). 

The sine-cosine transform enables us to delete redundant components from 
the real-valued data at negative frequencies. 


Consider a circular symmetric matrix, Ts , defined as follows: 


Ts 


c(0) c(l) c(2) c(l) 

c(l) c(0) ...c(3) c(2) 

c(l) c(2) c(l) c(0) 


(4.97) 


The circulant matrix Ts given above can be derived from Tc defined by eqn. 

(4. 44) by substituting C (N-1) = C (1), C (N-2) = C (2), The characteristic 

roots of Ts are (Fuller, 1976, pp. 135-137) 


(N-l)/3 

c(n) cos 2-rnj/N 

N /3 

Yj c ( n) cos 2TTn j/N 


for N odd 
j =0,1,...,(N-1) 
for N even 


(4.98) 


The equation above has a root for j = 0 and (N-l)/2 repeated roots for 
j = 1, 2, . . . , (N-1)/2, and finally another root for j = N/2 when N is even. We 
can find two orthogonal vectors for each of these repeated roots. These vectors are 
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11 


!»«••# 1 1 for j “ 0 

aj = ^[i,cos2it^,cos2tt^ cos2n^^^^], j = 1,2 NH(4.99.a) 

aj = ^^[1, -i, 1,-1,...,-1] j = N/2 and N even 

bj = [o, sin 2 n^, sln2n-^, , . . , sin 2 tt , j = 1,2, . . . , Jsh (4.99.b) 


where 


N„ 


r(N-l)/2 for N odd 


t N/2 - 1 for N even 
Hence we can form an orthogonal matrix Q as given below 

/S’ rA' 


where 


^ /2 rA"i 

^ = 7F LbJ 

( r"aQ//S 


(4. 100) 


A=< 




aw„ 


if N odd 


ao/ /S’ 


a. 


and 


B = 


aw/a 

bi 

b 

m 

tVi,, 


If N even 


Since is composed of the N characteristic vectors given by eqns. (4. 99. a) 
and (4.99.b) we can write 


Q Ts Q ~ Diag ( A.q , , . . . , ) 

where Xj is defined by eqn. (4.98). 


(4. 101) 


There exists a linear relation between the characteristic vectors of Q and 
the characteristic vectors of the F matrix defined previously. The F-matrlx 
given by eqn. (4. 30) can be written as follows 


fo 

F = I 

L 

where fj = ^ [ 1, , . . . , ] 
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From eqns. (4.99.a-b) and (4. 101) it Is possible to write 
1 

, llj IN-JJ 

j = 1,2,. ...N, 


aj = 

bj = 


75 [fj + 
[fj - 


Using the relation given by eqn. (4.102), the Q matrix is written 
Q = H F 


(N-l)/2 (N-l)/2 


where 


H =< 


f 

&! J 

1 

7? 

I 

.1 


i‘--. 

7l 

• 

< 

_ i ‘‘-i 

^ 1 

• 

-i-‘ 

1 

-1 

1 

1 1- 1 

I .-1 

/? 

1 • . 1 
1 *1 

1 

! 


___! ■•il_ 

Ir 

1 

! ! y^\ 


1 i 

1 •. I 

I Ti 

1 

1 

L 1 _ '-ii 

|-l*_ ^ J 


(N/2)-l 

(N/2)-l 


(N-l)/2 

(N-l)/2 

(N/2)-l 

(N/2)-l 


(4. 102) 


(4. 103) 


N = odd 


(4. 104. a) 


N = even 


(4.104.b) 


Now we can write down the frequency domain equations (windowed) when 
the sine-cosine transform (Q) instead of the Fourier transform ( F) Is used 
for the computation of y ■” Cxx X . The x-data vector is transformed into the 
frequency domain by 


X = Ag X (4.105) 

and the y-solution vector (unknown) by 

Y = (At)‘^ y (4.106) 

where Ag is the transform matrix defined as 

Ag = H F w ,4 107) 

For the covariance matrix of X given by eqn. (4. 105) we have (analogous to the 
derivation of eqn. (4. 35)) 

Cxx ~ Ag Cxx A^ 

= H F w Cxx W f'*’ H+' (4. 108) 
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or equivalently 

Cxx = 2 K HaN Fsn w 6xx w Fat hJ (4. 109. a) 

where K, FgN» w, and txx are defined by eqns. (4. 77.a-d) respectively, and 
Is defined by eqn. (4. 104,b) after letting N -* 2N. Eqn. (4. 109. a) can also be 
written as 

Cxx ^ 2 K Hsn W Fsn Fsn Cxx Fsn Fsn W FgN FgM K 
= 2 K HaN ^ Cxx K ^ 

A A ^ A A 

Wh©re W = FgN W FgN » Cxx “ FgN Cxx FgN 

If we denote 

Wh »= HgN ^ (4.110) 

then we can write 

Cxx = 2 K Wh Cxx wJ" (4. 111. a) 


(4.109.b) 
(4. 109. c) 


Eqn. (4. 111. a) is slmilaij^to eqn. (4.78), the only difference being the replacement 
of W of eqn. (4. 78) by Wh in eqn. (4.111. a). Although ^ of eqn. (4.78) is a 
circular matrix, Wh of eqn. (4. 111. a) is not. Therefore, we will use eqn. (4. 109.b) 
to determine the elements of Cxx • Through the substitution of the definitions for 
K, Hn , and Han it is easy to see that 


K Hsn = K 

Thus we can write eqn. (4. 109. b) as follows 

Cxx = 2 Hn K W Cxx W (4. 112) 


Now substitute jpqns. (^. 104. b) and (4. 77. a) for Hn and K respectively and eqn. 
(4. 109. c) for W and Cxx in eqn. (4.112). After some simplifications, we obtain 


SN-l 


^ for 0 ^ j,k ^ N/2 

^U0e-3J + %-3N+3j)(Pt-«3k -3T ^jt+N-3k)^i lN/2<k^(N-l) 


_ - 7o“'.V^N/2' (4.113) 

^7^^jt*-N-3j “ -3N )(^jt-3k'*' lN/2<j S(N-l) 

Xj(^f.+N-2J ~ -3N )(^^ + N-3 ir^X,+ 3k-3N N/2 < j , k ^ { N-1) 

V jt=o 


where n„ and co. are given by eqns. (4. 81. b) and (4. 82. b) respectively and 
f2_n = fiskN-n * K = 0,l,2, .... Since Cxx given by eqn. (4. 113) is approximately 
band-diagonal, as shown in Figure 4.7 (Heller et. al. , 1977, p. 41), we only 
compute the elements of the band-diagonal, so we can write 
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-3N49j)(^^-a(j+ ■) + Pt-aN-ts (j+i) ) 

for 0 s j s n/2, 0 s m Nb such that m+j ^ N/2 

(4. 114) 

3N-1 _ 

^E^9tr»-N-3j“ 94«J -*N)(^j^-N-a(j+, ) “ 9 ^(J+b)-3N ^ 

for N/2 <js:(N-l), 0 ^ m s^Nb such that m+j « (N-1) 
where Nb is the width of the baud. 

Nb + 1 



Figure 4. 7. Structure of Cxx under sine-cosine 
transformation (N = even). 


Having computed Cxx we can go back to the solution of the system of linear 
equations 

Cxx Y = X (4.115) 

where X, Y, and Cxx are given in eqns. (4.105), (4.106), and (4. 114), respectively. 
The equation above can be solved for Y by using the standard band-diagonal Cholesky 
decomposition which numerical computations are proportional to (Nb^ N). Finally, 
we can compute the space domain solution vector (y) from eqn. (4. 106) to obtain 

y = C7x^ X = w F+H'^’Y (4.116) 

Thus, the windowed frequency domain solution using sine-cosine coefficients, in 
order to compute the s-signal vector, is completed. 

Since we use a highly tapered window, such as the Kaiser window, in order 
to minimize the width of the band, Cxx will have very small eigenvalues. In 
other words, Cxx is almost ill-conditioned (Heller et. al. , 1977, p. 42). Con- 
sequently, Cxx has to be modified before solving eqn. (4. 115). This modification 
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can be done by adding a small amount, say 5, to the diagonal elements of Cxx 
to obtain 

Cxx = Cxx + 6l (4.117) 

Thus, Cxx cannot have any eigenvalues smaller than 6, therefore, Cxx is 
well-c ondltioned. 

The addition of (6 I) to Cxx actually corresponds to the addition of 

D = 6 • Dlag ( w?, w^, . . . , \v ^~^ ) (4, 118) 

where Wj are the window coefficients, to the space domain covariance matrix, 

Cxxf (Ibid., p. 43). A detailed study on the introduction of 6 and its consequences 
can be found in Heller et. al. (1977, pp. 42-47). 


4.6 Frequency Domain Collocation on a Two-Dimensional Region 
4. 6. 1 Large N and Two-Dimensional Wiener Filtering 

In previous sections, we have considered equally spaced observations along 
a single profile exclusively for the determination of the s-signal vector. In this 
section, we will consider observations given along more than one profile leading 
us to the two-dimensional frequency domain collocation. In order to obtain 
covariance matrices (for observations) of block-Toeplitz form, we will select 
the profiles in such a way that observations form a grid as well as they are 
equally spaced along every individual profile. 

Let us denote that: 


M is the number of parallel profiles with N observations along each 
Xd is the data vector of (M N) observations 

CxpXp is the covariance matrix of observations (it has block-Toeplitz form) 
So is the signal vector desired 

CspXp is the cross-covariance matrix between signals and observations. 


We can transform xb into the frequency domain as follows 


where 


Fo = 


Fo 

Xo 


(4. 119) 





F 

0 .. . 

.. 0 


0 

• 

• 

• 

F... 

.. 0 

(4. 120) 

_0 

0 .. . 

.. F_ 

[(MN) X (MN)1 


here F is (NxN) Fourier matrix given by eqn. (4.30). 
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Hence the covariance matrix of X0 can be expressed as (analogous to the 
derivation of eqn. (4. 35)) 


CxoXp - Fo Cxj,Xq Pt 

Since Cx X has a block -Toeplitz form, we can write 

D D 


^ D*0 


Tn Tia 
Tai %0 Taa 

Tmi Th3 Tm3 


Txm 
la M 

Tw 


(4.121) 


(4. 122) 


where Tjj 
becomes 


are (N x N) square blocks of Toeplitz form so that eqn. (4. 121) 

FTnF;f FTiaF^... FT^, 

FTaiF^ FTaaF"^.... FTaHF'*' 

FTmxF'*’ FT«aF'’’. ... FThmF"^ 


(4.123) 


We demonstrated before that the elements of (F Tn F ) converge to the elements 
of the diagonal matrix 


Dll = 2 tt Diag [ tii(o{]), til (cox), . . . , til (C0 n-i)] (4.124) 


where tn(Wn) is the spectral density of Tn, and it is given by eqn. (4.47) evaluated 
at n = 0, 1, . . . ,N-1. 

Fuller (1974, pp. 308-310) shows that the elements of (F Tjk F^, j 7^ k , 
also converge to the elements of the diagonal matrix 

Djic = 2 tt Diag[tjk(cOo), tjic(COi), ..., tjic(CON-i)]. j k (4.125) 


where tjk(cOn) is the spectral density of Tjn , and is given by eqn. (4. 47) evaluated 
at n = 0, 1, . . . , N-1. Thus, we can write eqn. (4. 120) approximately as 


Cx^Xp- 


Pii Dia Dim 

D^ai Daa Daw 

Dll Dw 


(4.126) 


where Djk, j,k = 1, 2, ..., M, which are square sub- matrices of Cx^x t are 
diagonal. 


The solution of the Sq - signal vector from two dimensional data can be 
written (in the space domain) from eqn. (4.22) 

C7\ Xb 

D 0 


Sd - C* X 
0 0 
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(4. 127) 


In the detorminatlon of So above, the costly and time consuming part Is the 
Inversion of or equivalently the computation of (Cx“\p«D)- Let us denote 

yo = Cx" xp Jib (4. 128) 

and define the frequency domain representation of yo (c.f. to eqn. (4. 119)) as 


Yo = Fo yo (4.129) 

then we can write 

Yo = Cx"\ Xj (4.130) 

D D 

Since Cx(,Xo defined by eqn, (4. 126) Is formed by (N x N) diagonal sub-matrices, 
we can easily and cheaply Invert and compute the product (crj x Xp ). Then, by an 
Inverse Fourier transform, we obtain 

yo = Fo"*" Yo (4. 131) 


Finally, we compute the So -signal vector 

Sb — Cs X yo {‘it 132) 

0 0 

The algorithm described above is equivalent to the "two-dimensional Wiener filtering" 


4.6.2 Small N and Two-Dimensional Windowed Frequency Domain 
Collocation with Sine-Cosine Coefficients 


If N, the number of observations along a single profile, is small then the 
(F Tjit F"^) submatrices of eqn. (4. 123) do not converge to a diagonal matrix. 
Therefore, in the case of small N, we will use the windowed frequency domain 
collocation in a two dimensional region as explained below, in order to obtain 
real -valued and band-diagonal submatrices for Cx^xp the orthogonal matrix Q 
of sine-cosine coefficients given by eqn. (4.103) will be used instead of the imagi- 
nary matrix F given by eqn. (4.30) in the transformation from the space domain 
into the frequency domain. 

Consider the definitions M, Xb, Cxjjxp» So , Cs x » as given in the pro- 
ceeding section. Here the transformations are applieci as follows. The x^-data 
vector is transformed into the frequency domain by 


Xd = Qd Wd Xb 


(4.133) 


126 


where 


Q) = H) Pb = 


HP 0 0 

0 HP 0 


HP 


(4. 134. a) 


[(MN) X <MN)1 


Wo = 


w 0. 
0 w. 

0 0 . 


0 

0 

w 


(4.134.b) 


((MN) X (MN)] 


here w is the window matrix (diagonal with diagonal elements 
(WotWj , . . . ,W(g_i ) defined by the window function). 


and yo -solution vector is transformed by 
Yo = [(Qo wo)^]"^ yo 

and finally, the covariance matrix Cxq by 

CxpXo = Wo Cx„x„ Wo^ Qo"^ 


(4.135) 




(4. 136) 


Substitute eqn. (4. 134. a) for Q) , eqn. (4. 134.b) for Wo, and eqn. (4.122) for 


Cxq xq to obtain 


CxnX 


D'O 


HPwTiiWpV HPwTiaWpV HPwTi„wpV 

HPwTaiWP'’’H‘’' HPwTaaWP%‘^ HFwTaMwP%'^ 


H PwTm 1 w p'^h'*' HFwTm aWP'^H'^ H FwTm m w p'*h'^ 


We can write for each square (N x N) sub-matrbc of Cx^x^ 

Bjk =HPwTjkWP%+= 2KH^NIiNWTJkWpJH^NK^ 1 s j,k s m 


(4.137) 


(4.138) 


where K, H^n • F 3N t w are defined as in section 4.4, and 

is a circular matrix (2N x 2N) extended from Tjk as explained in 
section 4. 4. 2. 

Eqn. (4. 138) is band-diagonal as shown in section 4.5. The elements of the (N x N) 
band-diagonal matrix Bjit can be computed as in eqn. (4.116). 

As before, we are here interested in solving the system of linear equations 

Cx^xp Yo = Xo (4.139) 

Since Cx^x^ consists of sub-matrices v;hich are band-diagonal, eqn. (4. 139) can 
be solved by a modified Cholesky’s decomposition which requires (M® Nb^ N) 


- 127 - 


operations comparcci to (M'^ KT*) In the space domain solution for yo. Following 
the computation of Yo from eqn. (4, 139)we can find the yo-aolutlon vector from 
eqn. (4.135) 

yo = Wo Q)'*’ Yo (4. 140) 

Finally, we can wrtie for the Sb-signal vector 

So ^ Csp xq Jt) (4. 141) 

Thus, the windowed frequency domain collcxiatlon in two-dimensional region is 
completed. Due to extensive programming needed for this technique we did not 
apply this algorithm. 


4.7. Summary 


In this chapter we have covered one- and two-dimensional frequency domain 
least-squares collocation. For a fast solution we have introduced an approximation 
and neglected off-band diagonal (off -diagonal in the case of Wiener filtering) terms 
of the frequency domain representation Cxx of the covariance matrix, Cxx, of 
observations. These neglected terms cause an error in the prediction of signals. 

In order to minimize this error, and yet preserve the speed of the algorithm, 
effective wb are used so that the off -band diagonal terms of Cxx » the frequency 


domain cc /. 
an optimup’ 
imatiOT nei 


natrix of tW windowed data, are as small as possible. Designing 
or a rarticular purpose and the error analysis due to the approx- 
oove are open for future studies. 
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< onclufllonu 


As we stated before, the modem Instrumentation cnabl(*d us to eolIe(!t a 
large amount of geodetic data. In order to process and Interpmt this dato 
effickuitly we Introduced the frequency metho<l8 and the fast Inversion of 
Teoplitz matrices. 

In (’hapter 2 we reduced the complicated function of geoid heights along a 
single profile to a series of simple trigonometric functions in order to Investi- 
gate the resolution (or cut-off frequency) of the (iK()S-3 altimeter data us a 
function of wavelengths or frequencies. The minimum full wavelength recov(^ra’ole 
from CjKOS- 3 altimeter data has been found to be about 100 km with the assumption 
of 27 mgals standard deviation on predicted ix)int gravity anomalies. We also 
computed the t.jtal power in the sea surface topography at and above 20 cycles/ 
revolution, which is equivalent to the total power with respect to the (JEM-O 
surface. This 'X)wer is about 2. 1 meters from some 9 profiles described in 
Table 2,3. 

In Chapter 3 we introduced a rigorous fast inversion algorithm for matrices 
of simple or block Toeplitz forms. In least-squares collocation solutions we have 
to invert the covariance matrix of observations (denoted by C = C*x throughout 
the paper). This covariance matrix has a dimension, say N, equal to the number 
of observations. In the case of equally weighted and spaced observations along 
a single profile and an isotropic and global (stationary) covariance function used in 
the computation of C , the resulting covariance matrix, C , is of simple Toeplitz 
form. The Inversion of such a covariance matrix requires a numerical operation 
proportional to N* compared to N® in case of conventional inversions. If we 
have observations forming a Cartesian grid with the properties explained above, 
then C is of block- Toeplitz form. The inversion of such a block Toeplitz matrix 
requires rou^ly (N® P®) numerical operations, where N in equal to the number 
of observations along a single profile and P is equal to the number of profiles, 
compared to (N^ P^) in the case of classical inversions. Another advantage of 
Toeplitz inversion algorithms over classical inversion algorithms are that we only 
Iiave to store one row in the case of simple Toeplitz matrices and one-block row in 
the case of block Toeplitz matrices. The only disadvantages in the application of 
Toeplitz algorithms are the following; 

a) The requirement of gridded data consisting of equally weighted and spaced 
observations along, at least, a certain direction. 

b) The requirement of an isotropic and global (stationary) covariance function 
In the computation of C. This Is also assumed In regular collocation. 

In our computations of gravity anomalies from geoid heights we used the theoretical 
covariance function given by the subroutine COVA of Tscheming and Rapp (1974). 
This function satisfies the second requirement, to order to satisfy the first 
requirement we created gridded data from GEOS-3 altimeter data using the 
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subroutine* PREDICT of the* Geodetic Soie>nce Plotting Package (GSPP) written by 
Hiinkel (1979). By using the subroutine PREDICT, It took roughly 1 second CPU 
time on Amdahl 470 of (XSU for the prediction of 1000 points. The mean variance 
of original observations was assigned to each of the predicted data. These are 
the ouly approximations introduced when we compare Toeplii'^. solutions and 
rigorous least-squares solutions. So it is up to the user to decide between the 
gain in computer time, storage and approximations. In our cases, as we showed 
in Chapter 3, the ai^roxlmations are negligible. 

The simple and block Toeplitz algorithms are much more efficient than the 
classical inversion algorithms, but for very large N, say greater than 1000, 
they are time consuming and become inefficient as well. Yet another faster 
method is the frequency domain collocation discussed In Chapter 4. Here every 
quantity in the space domain is transformed into the frequency domain and the 
solution in performed there. For a fast solution we take advantage of the simple 
structure of the frequency domain covariance and cross-covariance matrices. In 
the limit as N ® these matrices become diagonal and the method reduces to 
the generalized Wiener filtering. By the use of the frequency domain least-squares 
collocation, a considerable gain in computer time and storage Is obtained in com- 
parison with conventional least-squares collocation or least-squares collocation 
with Toeplitz inversion. However this method has disadvantages as well, such 
as the aforementioned disadvantages in the application of Toeplitz inversion al- 
gorithms and the errors due to edge effects. 

The methods presented in Chapters 3 and 4 shov> that they are effective tools 
for a fast calculation of desired signals from large amounts of gravimetric data, hi 
previous chapters these methods have been demonstrated for obtaining gravity 
anomalies from geoid heights. 
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Appendix l.A 


Fourier Series Representation of Periodic Functions and Parseval*s Theorem 


A periodic function x(t) with period T, i.e. x(t)=x(t + T) can be 
arbitrarily closely approximated by a harmonic polynomial of degree N and 
period T such that 

N 

(^n cos nWot + b„ sin ncOot) (l.A.l) 

The closeness of JiN(t) to x(t) is usually computed by least-squares approxi- 
mation over the interval ( 0, T ) defined by 

v^(t) = Jjx(t) -XN(t)]®dt (l.A. 2) 

V® ( t) , whose parameters are ag , % , b^ „ . . . , a(g , hj , is to be minimized. So 
let us Substitute eqn. (l.A. 1) in (l.A. 2) to obtain 

r ^ 

v®(t)= [x(t)--^- Y. (a„ cos nujot + bn sin ncogt)]® dt (l.A. 3) 

^ n=l 


Now let us take the dei’ivatives of eqn. (l.A. 3> with respect to the parameters 

= 2 f[x(t)--|s- (a„ cos ncoot +b„ sin nWot)]cos mofet dt (i.A.4.a) 

na, ^ nm\ 

= 2 I [x(t) --|2- y (a„cos ncogt +bn sinntOot)]sin mo^t dt (i.A.4.b) 

D, •» T Z 

The solutions to Sv®(t)/3a, = Sv®(t)/Sb, = 0 give the desired parameters. 

Recall the orthogonality relations defined by eqns. (1.2.a-b) in order to solve 
eqns. (l.A. 4. a-b) for the parameters 

3-n = "Ir f x( t) cos ncogt dt , n = 0,l,...,N (l.A. 5. a) 

T Jt 

bn = -^ J^x(t) sinnojotdt n = l,2,...,N (l.A.b.b) 

How close does eqn. (l.A. 1) ajproximate the true function x(t) ? hi order to 
answer this question open the parenthesis of eqn. (l.A. 2) to get 

v^(t) = x®(t) dt + Xi^®(t) dt - 2j^ x(t) JiN(t) dt (l.A. 6) 

Using the orthogonality relations mentioned above we can easily derive thiti 

Y «N®(t) dt=-|r|^ [-|“+ (a„cosncOot+bnSinno3ot)]^dt= (yj + i (an®+bn®) 

(l.A. 7) 

(Eqn. (l.A. 7) is called ParsevaUs Theorem ) and 
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HHtCEuiuu hMiac olmihK NOT FILMED 


(1.A.8) 


Therefore eqn. (1. A. 6) can be written as 

r 8 T> ^ 

v^(t) = dt - lTy-^)+ + a 0 (1.A.9) 

or 

” ^ [(^) + + (I.A.IO) 

Eqn. (1. A. 10) is known as Bessel's inequality. From this Inequality we can read 
that 

lim a„ = lim b„ = 0 

0 -^00 Q -400 


Now let y(t) be equal to Xig(t) in the limit as N-*™, i.e. 

00 

y(t) = lim XN(t) = ^ + E (a„cosnatet + b„3innoJbt) (l.A.ll) 

where a„ and b„ are Fourier coefficients of square integrable function x(t). 

By Bessel's inequality the partial sum of squares of the Fourier coefficients 

l^ryw* nwyxr ormn^^a • ^ 4- J rnt«^ 4.1. iil.. .4.2— « 4. \ 

3'-’ Q j.vy* uiJiy £>v|u.a.4.v/ iuiAJgXixiL/xo xuu«^uxv^ii« X xic; j. C7 J.U Tc; uuc xuuCtlOu y(l) is 
square integrable on (0,T), Now consider 


d(t) = x(t) -y(t) 


(l.A. 12) 


then, d(t) is also square integrable on (0,T). By Shwarz' inequality we can write 
write 


0 s Jjd(t)|®dt s J^|d(t)Idt (l.A. 13) 

This expression is equal to zero by the "theorem of uniqueness of integrable 
functions". That is to say 


d(t) = 0 , x(t) = y(t) 


(l.A. 14) 


Thus all the Fourier coefficients in the Fourier expansion of d(t) are zero 
according to the "theorem of uniqueness of continuous functions". Hence we can 
conclude that 

an °° 

x(t) = limxg(t) “ 2 ^ ^ ( a„ cos nofet + bnSln ncOot) (l.A. 15) 
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Appendix l.B 

Fejer Kernel and the Derivation 


Substitute eqns. (1.3.a-b) in (1.44) to obtain 


jf{ J]^l-^;^^(X)lco8 note X cos nccbt + sinntObXsinncOot] 

j-dX + I^J x(X) dX 


(l.B.l) 

Now let z = ote( X - 1) to have 



(l.B. 2) 

Mt) = "t Jt x(X)dX 

(l.B. 3) 


. N 

where Fn+i(z) = 1 + 2^2i(l -u/(N+l) ))*cos nz is called the Fejer kernel. The 
Fejer kernel consists of cosine functions as we can see from the definition above. 
Thus it is possible to express the Fejer kernel as a sum of cosine functions as 
follows: 

Fn+i(z) = 1 + 2^1 -^j^osz + 2(l-~^os2z+ ... + 2^1-jj^cosNz (l.B. 4) 
Multiplying eqn. (l.B. 4) by sin z/2 from both sides we obtain, 
sin|- FN+i(z) = sin|- + 2 (l-^^^sin|-cos z + ... + 2^1-^^:^8in|- cos Nz (l.B. 5) 
Using the trigonometric identity 


2 sinx cosy = sin(y + x) - sin (y-x) 
we can easily prove that 

sin| F«.,( z ) = sln| + (l-ji)(sln|z -8taiz) + (l- + . . . 

■(l-^^)(sin(N+^)z-sin(N-i)zj 


(l.B.b) 


1 • z 1 3 N . 


(l.B. 7) 


Thus the Fejer kernel can be expressed as 

N 

Fn+i(z) = X! siin ^ I 


1 ^ sin ( n z ) 

N+1 „eo 


(l.B. 8) 


The expression for the Fejer kernel is not used in its present form. It can be 
simplified further, first by writing 
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(1.B.9) 


X X ' z 

“ iw iFTT.So “‘”1 “ 

then, opening the sum of the expression above to get 

Fn+i(z) = [2 sin^z sin^z + 2slniz sin|’z+. ..+2slnizsln(N+i)z] 

finaUy using the trigonometric identity 

2 sin X sin y = cos(x - y) - cos(x+y) 


we obtain 


Fn+i( z ) 


_ sin ^ jz 

2(N+1) 
_ sin~^ jz 
2(N+1) 
siif^ jz 
2 ( N+1 ) 


[2-2 cos®i z + cos z ~ cos 2z + cos 2z + . , , + cos 
[2-2 cos° ^ + cos z - cos ( N+1 )z ] 

a 

[ 1 - cos® (N+1 )| + sin® ( N+1 ) I ] 

2 2 


(N 


The equation above is equivalent to 


(I.B.IO) 


(l.B.ll) 


l)z-cos(N+l)z ] 


(1.B.12) 


p = JL r-S in(N+l )L z/2)-] ® 

Fn+i(z) gin(z/2) J 


(1.B.13) 


Appendix l.C 

VarianoeB of Fourier Coefflolenta 


Let us suppose that the function x(k*At) Is Independent and normally 
distributed N(/i , and has the Fourier coefficients defined by eqn. (1.93), 

? x(k*At){ , n=0,l,...,NH (l.C.l) 

\hj N uko I sln2nkn/NJ t » » h v 

M - / (N/2) If N even 
where N„ - | ^ 

So If we denote 

Y = [ao, aj X = [x(0).x(At) x((N-l)At)]^ 


=N 


1 

2nl 


1 ... 

4tt1 


1 cos cos .. 

N N 


.cos 


2(N-l)TTi 


N 


1 ooeiiia „03 2('n-1)wNh 

N N N 


then we can write 


(l.C. 2; 


Y = G X 


(l.C. 3) 


The variance -covariance matrix of the function Y can be e^ressed as 
follows (UotUa, 1967): 

E = G g' (l.C. 4) 

VJy X_» jj 

where Ev Is the variance-covariance matrix of the Fourier coefficients (Uq, 

51.. ..,aN„) and Ex is the variance -covariance matrix of x(k*At), k=0, 

1. . . . , ( N-1 ) , which Is diagonal with diagonal elements equal to by assump- 
tion, l.e. 


- diag (o^, a® o^) 

Therefore, we have. 




n T 

CT^ G G 


(l.C. 5) 


(l.C. 6) 
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We know already that G and are orthogonal such that 


G G^ 


4 


N 0 0 


2 0 0 

0 (N/2) 0 

. J. 

0 1 0 

• • • 

• • • 

• • • 

• • • 

N 

• • • 

• • • 

0 0 (N/2) 


6 0 i 


Thus eqn. (1.C.4) can be written as 


cr?„ = var(a„) 



for n = 1, 2, . . . , (Nh -1) 

for n = 0, or n = Nh and N even 

for n = Nh and N even 


Similarly we can prove that 


<^b„= var(b„) 



for n = 1, 2 , . . . ,(Nh -1) 
for n = Nh and N even 
for n = Nh and N odd 


(1.C.7) 


(1.C.8) 


(1.C.9) 


- 140 - 


Appendix l.D 

The Fast Fourier Transform ( FFT) 


The number of operations required by the conventional Fourier transform 
method is proportional to , N being the number of data points. On the other 
hand, the FFT algorithm requires a number of operations, which is proportional 
to (N loga N) . The FFT is an algorithm by which the discrete Fourier trans- 
form (DFT) can be computed much more rapidly than by other available algorithms. 

Let DFT be defined as in eqns. (1. and (1. 66), namely, 

"E x(k-At) (l.D.l) 

N k aO 

x(k-At) = V X„ (I.D.2) 

DM 0 

When N is a product, say N - r*s , then the Fourier transform can be calculated 
in a two-stage process. That is to say, as if Xn and x(k*At) were defined on 
two dimensional (rxs) arrays with array indices (k^ , ko) and (Ui , no) such that 

k > (ki,ko) (l.D. 3) 

k “ kx • r + ko for ko ~ 0, 1, . . . ,(r-l) and kx = 0,1, . . , ,(s-l) 

n > (nx,no) (l.D. 4) 

n = nx* s + no for no = 0 , 1 , ... , (s-1) and nx - 0, 1 , . . . , (r-1) 

we define Wn = exp(2rri/N) (l.D. 5) 


and use e<ps. (1. D. 3) and (1. D. 4) , then we can write 
Wn’'" = Wn‘‘i''i”« Wn'*''”!* • Wn“i“o' • 




(1. D. 6) 

Thus e(]pi. (l.D. 2) becomes 



• -1 r-1 


x(k» At) 

= x[(kx,ko)-At] = E E X(„ . )-w,‘‘°“i.w,''i”o-w,'‘<=Po 

ngSsO n^sO ' 

(l.D. 7) 

If we denote 

r-1 



iw H o 

(1. D.8) 


then, 

• -1 

x[(kx,ko)*At] = ^ico,no)' 
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(1.D.9) 


Notice that first we sum over % in eqn. (1. D. 8) to form an intermediate 
array ^^^Uo*nQ) and then sum over no . Hence X\ko»«o) Wn‘‘o times 
a set of r-tenn Fourier transforms. Finally eqn. (1.D.9) Is a set of s-term 
Fourier series with as coefficients. The total computation amounts 

r». s + s® • r = N ( r + B ) (1. D. 10) 

for this particular two-stage process. 

Similarly, it can be shown that it N =» r^ x r^ x . . , x r, , then the number 
of operations is equal to 

N X ( rj + rg +. . . + r, ) 

If we further assume that rj = Tg = ... = r,, then N = r" and the number of 
operations is equal to 

N»m«r = N(log,N)r = N(log^ N) r/logj,r (l.D.ll) 

when r is equal to three, we obtain the lowest number of calculations. This can 
be proven by minimizing the above equation, i. e. 

y ( r ) = N (logg N) r/logg r to be minimized (1. d. 12) 

The solution to y ’( r ) = 0 gives the optimum factor r . 

9 

y’(r) = N(loggN) [loggr - r (logg r) ] /(logg r)® » o (1.D.13) 

This is equivalent to 

a 1 

logaT- r^(loggr) = log^r--^ = 0 
The solution of the above equation iS; 

r = 2 • « 3 (nearest integer) (1, d. 14) 

However for r = 2, 4, 8, . . . we can further reduce calculations by avoiding 
multiplications when the powers of w^ are simple numbers like ±1, ±1. There- 
fore the most commonly used FFT algorithms use r = 2, i.e. N = 2" . Generally 
a seperate derivation is given for this case. Let 


x(k*At) 

x„ 

XJ I 
x:./ 

k,n = 0,1,..., (N-l) 

(1. D. 15) 

x[(2k'+l)* At] > 

k'n' = 0,l,...,((N/2)-l) 

(1. D. 16) 


then the discrete Fourier transform (DFT) can be written as 

j N-l 

Xni = “ ^ x(k»At) Wn“’"‘’ (l.D. 17) 

^ icttO 

using eqn. (l.D. 16) we have 
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(1.D.18) 


Xrt’ = (x(2k'* At)*w,7"'‘''‘' + x((2k'+l)*At)‘WN"^*'‘'*^ 

N ^ k* a«0 

But we know that w,® " , therefore 

(N^>-l 1 (n/^- 1 -k'n* -nS 

= *W) Xx(2k’*^t).w,^ ^t(2k'+l)-Atl-WH/3-w.| 


which can be written as 


and 


x„' = i [ xj + xj' Wn"’'] 

X„»+(N/a) •-= i IX„'« - Xrf Wn""'] 


(1.1) . 19) 

(1. 1) . 20) 
(1.D.21) 


Notice that X’nt and X'*„i in eqn. (l.D. 21) ar«j also calculated through 
’’doubling". So this successive doubling continues until it is no longer a multiple 
integer of 2 . 

When the FFT algorithm is used, we avoid comiputations by a factor pro- 
portional to cKm , where 0 (n is u constant times the ratio of the number of opera- 
tions in the conventional method over that of FFT, i»e. 

a„ = c N®/(NlogaN) = cN/logfcN (1. D. 22) 

where c is a constant greater than one . The computational reduction by a factor 
0 !n (assuming c = 1) for various N are tabulated below: 


N «N 


2 

2 

16 

4 

256 

32 

1024 

341 

4096 

341 

8192 

630 

16384 

1170 


Since less operations are required in FFT computations, the truncation 
error is considerably smaller than that of conventional computations. Thus FFT 
gives more accurate results. 

In order to demonstrate a simple application of the FFT algorithm, the 
discrete data given at N points were transfonned by the "FFTR" subroutine 
of the IMSL library and by the conventional algorithm. The results are given 
below: 
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N 


Conventional Method (sec.) FFT Method (sec.) 


256 2.52 0.06 

512 9.92 0.13 


The time required for computing Fourier transform by conventional al 
gorithms and by the FFT algorithm is lUustrated in Figure l.D. 1 below (ex- 
tracted from Cooley, 1969). 



Figure 1. D. 1. Time required for calculation of Fourier transform 
of real data on IBM 7094 using FORTRAN with con- 
ventional and fast methods. 

We have tried to introduce the FFT algorithm as a computational tool here. 
For detailed discussions see Brigham (1974) and Path (1974). 
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Appendix 3. A. A FORTRAN Program for the Inversion of Simple Toeplltz Matrices 


c 

HDBROUTIIfE TOEPL(B,C,X.F.MAX) 




A 

10 

c 

0 

c 





A 

A 

20 

FUNCTION - TO COMPUTi: THE TOEFL ITC 

INVERSE 

OF ^T’ 

A 

A 

40 

c 

AND MULTIPLY IT BY THE DATA VECTOR *X* 

A 

50 

c 

TO GET *F* 




A 

60 

c 

PAiL\METERS B - THE INTERMEDIATE VECTOR OF LENGTH -MAX- 

A 

70 

G 

TO GET TBE INVERSE 




A 

B0 

c 

C - VECTOR OF LENGTH -MAX- .CONTAINS 

THE FIRST 

A 

90 

G 

'T' ON INPUT AND REPLACED UY 

BY INVERSE 

A 

100 

C 

ROW OF »T’ ON INPUT AND 

REPLACED BY THE 

A 

1 10 

G 

INVERSE LATER 




A 

120 

G 

X - DATA VECTOR OF LENGTH -MAX- 



A 

130 

C 

F - VECIOR OF LENGTH -MAX- AND 13 

PHODUCl OF 

A 

140 

G 

TINV X 




A 

150 

G 

MAX - THE DIMENSION OF THE TOEPLITZ 

MATRIX ’T* 

A 

160 






A 

i 

t j 





1 « 

G 





A 

100 


SUBROUTINE TOEFL ( B.C, X, F. MAX) 




A 

190 


IMPLICIT REAL?#:8<A-H,0-Z) 




A 

200 


DIMENSION B( 1) , C( 1) , X(l), F< 1) 




A 

210 

G 

• • • 




A 

220 


B(MAX)=1.D0 




A 

230 

G 

.. COMPUTE ’AA* AND ’E’ VARIABLES OF THE ALGORITHM 



A 

240 


DO 13 K=1,MAX 




A 

250 


E=0.D0 




A 

260 


AA=0.D0 




A 

270 


DO 11 L«1,K 




A 

280 


KA^riAX-L+l 




A 

290 


AA=AA+B(NA) 




A 

300 


E=E+B( NA) K-L+2) 




A 

310 

n 

CONTINUE 




A 

320 


AA= 1 . DO/AA 




A 

330 


IF(K.Ea.MAX) GO TO 13 




A 

340 

c 

GOM.PIJTE ’B’ VECTOR USING REGUPjSIVE ALGORITHM 




A 

350 


EAA--Esf:AA 




A 

360 


NB=rLlX-K 




A 

370 


B( NB)=EAA 




A 

380 


IFCC.Ea. 1) GO TO 13 




A 

390 


ICHALF=IC/2 




A 

400 


ICR=I'I0D(K,2) 




A 

410 


DO 12 LL-l.KHALF 




A 

420 


NB1=NB+LL 




A 

430 


riB2-MAl^LL 




A 

440 


T1 = B(NE1) 




A 

450 


B( NBl ) ^B( NBl ) +EAA*B( NB2) 




A 

460 


IF(KR.Ea.0.AND.LL.Ea.KHALF) GO TO 12 




A 

470 


B( NB2) ^B( NB2) +EAAJf:Tl 




A 

480 

12 

CONTINUE 




A 

490 

13 

CONTINUE 




A 

500 

c 

. . . COMPUTE THE LAST ROW OF THE INVERSE 




A 

510 


DO 14 1=1, MAX 




A 

520 

14 

C( I)=AA:#:B( I) 




A 

530 

G 

. . . MULTIPLY THE COMPUTED ELEMENTS OF THE INVERSE 

BY THE 

CORR.FART 

A 

540 

G 

OF X-VECTOR 




A 

550 


NI = MAX 




A 

560 


DO 16 NJ=1,MAX 




A 

570 

15 

GALL TMULT ( C, X, F, NI , NJ , MAX) 




A 

580 

G 

. . . COMPUTE THE REMAINING ELEMENTS OF THE INVERSE 

AND MULTIPLY BY 

A 

590 

G 

THE CORRESPONDING ELEMENTS OF THE X-VECTOR 




A 

600 


N=1 




A 

610 


NN=MAX 




A 

620 


M=MAX-1 




A 

630 


MHi\LF= MAX/2+1 




A 

640 


MM=n+NHALF 




A 

650 


DO 17 I = MHALF,M 




A 

660 


II = MM-I 




A 

670 


N=NM 




A 

680 


NN= NN- 1 ^ ; ; IV. ^ .... _ 

NL=I-MH4LF+1 , - r 

T2^AA^Bi lU 




A 

A 

690 

700 





A 

A 

710 

720 


DO 16 J=N,NN 




A 

730 


NBl=M-J+i 




A 

740 


T3=B( NBl):^cTl 




A 

760 


T4=B( J):fir2 




A 

760 


C( J)=C( J+D-T3+T4 




A 

770 

16 

CALL TMULT ( C, X, F, 1 1 , J, MAX) 




A 

780 

17 

CONTINUE 




A 

790 


RErORN 




A 

800 


END 




A 

810 




80BR0DTIWE TPTOLT (C,X,F,RI,IfJ»MAX) B 10 

. . . THIS SUBROUTINE CONSIDERS SVMHETRY, PERSYMMETRY AND COMPUTES B 30 

THE PRODUCT B 40 

...C(NJ) IS THE (NI,RJ)TH,<NJ,NI)TH,(NII,NJJ)TH,(NJJ,NII)TH ELE B 50 

OF THE INVERSE B 60 

NI = ROW NUMBER OF THE INVERSE B 70 


NJ = COLUMN NUMBER OF THE INVERSE 

IMPLICIT REAL?J^3(A-H,0-Z) 

DIMENSION CU; , X(l), F(l) 

NII = MAX-NJ+1 
NJJ=i..lX-‘NI+l 
F(NI)=F(NI)+C(NJ)*X(NJ) 
IF(NI.Ea.NJ.AlMD.NI.EO.WII) RETORN 
IF(NI.EQ.NJ) GO TO 1 
F(NJ)=F(NJ)+C(NJ):f:X(NI) 
IF(NI.Ea.NII) RETURN 
F( NJJ) =F( NJJ) *l G( NJ) W:X( NI I ) 

1 F(Nm=F(NII)+G(NJ)JftX(NJJ) 

RETURN 

END 


MATRIX B 60 

B 90 
B 100 
B 110 
B 120 
B 130 
B 140 
B 150 
B 160 
3 170 

B 160 
B 190 
B 200 
B 210 
B 220 


Appendix 3 .B. A FORTRAN Program for th^ Inversion of Block Toeplitz Matrices, 


c 

c 

c 

c 

c 

0 

c 

c 

c 

0 

c 

c 

c 

c 

V. 

c 

c 

c 

c 

c 

c 

0 

c 

0 

c 

c 


c 


8UBH0OTINE BTOFPU D, DUMMY, C, F., AA, EAA, T1 , TEMP, TO, BC,T», T4, AINV, IR, 
1C) 


FUNCTION 


PARAMETERS 


- TO OOMPirrE THE INVERSE OF BLOCK TOEPLITZ 
’T* AND MULTIPLY IT BV THE DATA VECTOR 
‘X* TO GET *F’ 

T - THE BLOCK TOEPLITZ MATRIX TO BE INVERTED 

H - THE INTERMEDIATE MATRIX OF DIHEN, ( IR, IC) 

TO OBTAIN Tin: INVERSE 

C - FIRST BLOCK ROW MATRIX OF ‘T* MATRIX ON 
LATER REPLACED BY THE BLOCK INVERSE ROWS 
INPUr,LATER REPLACED BY BLOCK INVERSE 
ROWS 

X - DATA VECTOR OF LENGTH -MAX- 

F - VECTOR OF LENGTH -MAX- .\ND IS PRODUCT OF 
TINV K 


JR - NO.O*' MAXIMUM ROWS POSSIBLE 
1C - NO. OF MAXIMUM COLUMNS POSSIBLE 
I MAX - DimiNSION OF THIC SQUARE SUBBLOCKS OF * T’ 
JMAX - DHIENSIOIn OF *T* ^WRIX 
MAX - RATIO OF JMAX OVER AX. MAX» JMAX/ IMAX 
( D A DUMMY ) CAN SILVRE THE SAME i.OCATIONS , . CALLING PROGRAM IS 
CAIX DTOEPL(B,B,C,i:,AA,EAA. Tl , TEM2\ TO . BG, TO , T4, AINV, IR, 1C) 


SimROUTINE BTOEPL I B, DU?IMY, C . E. AA. EAA, Tl . TO, BG, TO, T4 , AINV, IR, 

1 IC) 

ILLICIT ni5Al.?^:0( A-M,0-Z> 

DlPa:WSH)N B(IR, IC). DUmmiR, IC), C( IR. IC) 

DIMENSION E(IU, IR). AA( IR. IR) , EAA( IR. IR) , Ti( IR, IR) . TO( IR, IR) 
DIMENSION TO(IR, III), T4i IR. IR) . BG( IR, IR) , TEMPI IR, IR) , AINY(I) 
COM. ION . DlfffiN/ IMAK. JMAX, MAX 
l.OGICAI, TICK 

... INITiAIIZE -B- MATRIX TO ZERO 
CAU. D7A:U0 (B. IR5^JO) 

IMS - 

NBM-(MAK- I):MMAX 


A le 
A 20 
A 00 
A 40 
A 50 
A 00 
A 70 
A 00 
A 90 
A lOt* 

A no 

A 120 
A IBO 
A 140 
A 150 
A 160 
A 170 
A I no 
A 190 
A 200 
A 210 
A 220 
A 230 
A 240 
A 250 
A 260 
A 270 
A 2t>0 
A 290 
A 300 
A 310 
A 320 
A 330 
A 340 
A 350 
A 360 
A 370 
A 3«0 


11 

C 

4'' 


12 


13 

C 


C 

C 

C 


14 


DO n {-{.nvii: 

jj=rmM+i 

B( I,JJ)=^1.D0 


.. CONPirrE ’AA’ AND *E* MATRICES OF THE ALGORITHM 

DO 16 K=^i,MAX 

DO 12 1"J,IMAX 

DO 12 J--‘ I , I MAX 

E( l.J> =0.1)0 

AAl i.. 0=0. DO 

DO 13 L= :.K 

NBC= ( L- n I MAX 


NBA=(MAX-L)1cIMAX 

ID=2 

CALL WLT (D.C,TEMl’,TEIH%BG, IR, IC,NBA,NBC, ID) 

CALL ABS ( AA.BG.AA, IR, IC) 

IF(K.EO.MAX) GO TO 13 
NnC=(K-L+I)rn^lAX 

CALL mJLT (B.C,TEMP,TEm>,BG, IR, IC.NBA.NBC, ID) 

CALL ABS (E,BG.E. IR. IC) 

CONTINUE 

INVERT THE SQUARE PlATR IX ’iV\* 

CALL AAINV (AA.AINV.IR) 

IFIK.EO.MAX) GO TO 16 

... COWUTE ’D* MATRIX USING RECURSIVE ALGORITHM 
. . . *'rr IS THE THE SQUARE MATRIX TO BE ADDED TO THE PRODUCT 
’TO* TIMES ’EAA* .WIffiRF. EAA=E3^c AA 


ID= I 

CALi, ?niLT (DUiW.D^IMMY.E.AA.EAA, IR, IC. NX.NX, ID) 

pu:k=.falsk. 

DO 14 1 -..IMAH 
DO 14 J= I , I MAX 
EAA( I, J) “EAM I, J) 

NB=(MAK-K-1)>?UMAX 

CALL PUTIN (B.EAA, IR, IC.NB.PICK) 

IFIU.EO. 1) GO TO 16 ^ 

K11A!.F=K/:! 

PICK= .TIVTE. "V 

iai= M0D(K,2) 

DO 15 Url.KHALF 
NBI r-lH!.L:^:IMAK 

CALL PIIT.N (B.Tl, IR, IC.NBI.PICK) 

CALL rUT»N (R.TO, IH, IC.NB2.PICK) 


c 

V - ' 


A 390 
A 400 
A 410 
A 42J 
A 430 
A 440 
A 460 
A 460 
A 470 
\ 400 

A 490 
A 500 
A 510 
A 520 
A 630 
A 540 
A 560 
A 560 
A 570 
A 580 
A 69v9 
A 600 
A 610 
A 620 
A 630 
A 640 
A 650 
A 660 
A 670 
A 6B0 
A 690 
A 700 
A 710 
A 720 
A 730 
A 740 
A 750 
A 76f) 
A 770 
A 780 
A 790 
A 800 

A a 1 0 

A 820 
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c 

• • • 

A 

81^0 


1D» 1 

A 

840 


CALL MULT ( DUMMY, DUMMY, EAA,T2,T3, IR, IC.IlX.lfX, ID) 

A 

600 


CALL AABS (B,T1 ,T3, IK, IG.HFl) 

A 



IF(KR.Ea.0.AWD*LL.Ea.KHALF) GO TO 10 

A 

870 


CALL MULT (DUPIMY,DUMMY,EAA,T1 ,TT, IR, IC,KX,NX, 

A 

880 


CALL AABS ( B. *12,T3, IR, IC, NB2) 

A 


10 

CONTINUE 

A 


16 

CONTINUE 

A 

910 

C 

. . . COMPUTE THE LAST BLOCK ROW OF THE 

A 

920 


10=3 

A 

930 


BO 17 1=1, MAX 

A 

940 


NL=», I-1)J.^IMAX 

A 

950 


CALL MULT (B,C, AA, TEMP, TEMP, IR, IC,RL,NL> ID> 


960 

17 

CONTINUE 


970 


PICIC=.TRUE. 

A 

9G0 


DO 18 J=1,MAX 

A 

990 


NI=(MAX-l):faMAX 

A 

!000 


NJ^( J--1)><«IMAX 

A 

1010 


CALL PUTIN (C,T1, IH, IC,NJ,FICIO 

A 

1020 

IB 

CALL B^IULT (T1 , X, F, NI , NJ, IR. IC) 

A 

1030 

C 

* • > 

A 

1040 


N=1 

Ik 

1050 


NN=MAX 

A 

106> 


M=HAK“1 

A 

10V 0 


MHALF=MAX/2+l 

A 

1080 


MP=M+MHALF 

/ 

1090 


DO 20 I = l^mALF,M 

A 

1100 


II = MP--I 

A 

1110 


N=N-1 

A 

1120 


NN= NN-l 


1130 

C 

... COfiPUTE " T1 = B(N,N-S) TRANS? sic AA * 

A 

i lAi) 


ID=3 

A 

1100 


NL=(M-II>5?5IMAX 

A 

1160 


CALL MULT (B,DUMMY, AA,TEMP,T1 , IR, IC,NL,NL, ID) 

A 

117u 

G 

.. COWPUTF, ’ T2 = B(N, 8-1) TRANSPOSE ^ AA 

A 

1160 


NL=( II-l)sltIMAX 

A 

1190 


CALL MULT ( B, DUMMY, AA, TEMP, T2, IR, IC,NL,NL, ID) 

A 

1200 

G 

. . . COrIV’UTE THE REMAINING SUBBLOCKS OF THE JNV»^;E AND MULTIPLY BY 

A 

!M0 

C 

K VECTOR 

A 

1220 


10=4 

A 

1230 


DO 19 J= I,NN 

A 

1240 


NBl=(N-.T)JiiIKAX 

i* 



NB2=( J-1)::^IMAX 

A 



CALL MULT (B. DUMMY, T1 , TEMP , T3, IR, 1C, NBl ,NB1 , ID) 




CAL.-. MULT ( B, dummy, T 2. TEMP, T4, !R, rC,NB2.NB2, ID) 

A 



CALL ZZTi’ (C,T3,T4, IR, IC,NB2,TEMP) 

• 

*'^9o 

19 

GALL BMULT ( TEJIP, X, F, NL,RB2, IB, IC) 

A 


20 

CONTINUE 

A 

1 '>■» 


RETURN 

A 

1320 


END 

A 

:*^30 




is 
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„ .. OOOO DO ftOOO 


SUBROUTIHIS AAIJTV (AA^AINV* IR) 

. . . THIS SUBROUTINE COMPUTES THE INVERSE OF AA-MATRIX 
. . . AA IS THE MAIRIX TO BE INVERTED ON INPUT AND INVERSE ON OUTPUT 


IMPLICIT REAI.fO(A-H,0-Z) 
DIMENSION AA( IR, IR> , AINV( 1) 
COMTK)N /DIMEN/ IMAX, JMAX, MAX 
IN( I)^( LM-D/2 
N=^0 

DO I I«l, IMAX 
DO 1 J= 1 , I 
N=N+1 

1 AINV(N)=AA( I, J) 

CALL DSINV (AINV, IMAX, l.D-14,IER) 
!F( lER.Ea.O) GO TO 2 
WariE <6. 7) ISR 
STOP 10 

2 CONTINUE 

DO 6 1= 1, IMAX 
UO 0 J-tr IMAX 
IF(J-I) 3,4,4 

3 N=IN(I)+J 
GO TO 5 

4 N=IK(J) + I 

5 CONTINUE 

6 AA( I, J) = AINV(N> 

RETURN 


FOR^IAT ( lOX, ’ INVERSION ERROR lER » 
END 


B 10 
B 20 
B 30 
B 40 
B 50 
B 60 
B 70 
B 80 
L 90 
B 100 
B 110 
B 120 
B 130 
B 140 
B 150 
B 160 
B 170 
B 180 
B 190 
B 200 
B 210 
B 220 
B 230 
B 240 
B 25t:^ 
B 260 
B 270 
B 2B0 
B 290 
B 300 

14 ) B 310 

B 320 


SUBROUTINE PUTIN (A.B, IR, IC,NA,PICK) 


... -B- IS THE MATRIX TO BE INSERTED FOLLOWING (NA)TH COLUMN OF 
A-MAiaiX IF P ICK=. FALSE. , OR EI^E TO PICK FROM *A’ IF PIC&.TRUE. 

IMPLICIT REAL*t8(A-H,0-Z) 

DIMENSION A( IR, IC> , B( IR, IR) 

COMMON /OIMEIi/ N,M 
LOGICAL ?ICK 
IF(PICK) GO TO 2 
DO 1 1=1, N 
DO 1 J=1,N 
NN=NA+J 

1 A( I.NII)=B( I,J) 

RETURN 
DO 3 1=1, N 
DO 3 J=1,N 
NN=NA+J 
B( I, J)=A( I,NN) 

RETURN 
END 


G 10 
C 20 
C 30 
G 40 
C 50 
G 60 
C 70 
C 80 
C 90 
C 100 
C 110 
C 120 
C 130 
C 140 
C 150 
C 160 
C 170 
C 180 
C 190 
C 200 
C 210 


SUBROUTINE ABS C A,B.R, IR, IC) 


.. SUMS UP ^A’ AND ’B* MATRICES TO GET •R* MATRIX 

implicit REAL^8(A-H,0-Z) 

A(IR, IR), B(IR, IR), R( IR, IR) 
rf9i4._ 0IMEN/ N,M 
L*> . 

?>0 i 

i R( : I) 

rsTum 

END 


D 10 
D 20 
D 30 
D 40 
D 50 
D 60 
D 70 
D 80 
D 90 
D 100 
D 110 
D 120 
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SUBROUTINE AABS ( A»T,TT, IR» IG.NB) 

. . . SUMS UP T AND TT-MATRICES AND INBERTO IN A-HATRIX FOLLOWING 
( NB) TH COLUNBT 

IMPLICIT REALsf«8(A-H,0-Z) 

DIMENSION T(IR, IR), A( IR, IC) , TT( IR. IR) 

COMMON /DIMEN/ N.M.MAX 

DO* r I»1.N 
DO I J=l,N 
NN= NB+J 

A( i.NN)=T( i,j)+rr( I, J) 

RETURN 

END 


E le 
E 20 
E 30 
E 40 
E 50 
E 60 
E 70 
E 80 
E 90 
E 10U 
E 1 10 
E 120 
E 130 
E 140 
E 150 


SUBROUTINE MULT ( A,C,D,E,R, IR, IC,NJ1 ,NJ2, ID) 
IMPLICIT REAL:f:8(A-H.0-Z) 

DIMENSION A( IR, IG) , G( IR, IG) , D( IR, IR) 
COMMON /DIMEN'' N,M 


GO TO ( 1.3, 5, 7, 9) . ID 

CONTINUE 

DO 2 1=1. N 

DO 2 J=1.N 

R( I, J)=0.D0 

BO 2 K=l,N 

R( I. J)=R( I. J) fD( I.K):f:E(K.J) 

RETURN 

DO 4 1=1. N 

DO 4 J=1.N 

R( I. J) = O.D0 

NC=NJ2+J 

DO 4 K=:.N 

NA=NJ1+K 

R( I, J)=R( I. J)+A( I,NA)*C(K,NC) 

RETURN 

DO 6 1=1. N 

DO 6 J=l,N 

JJ=NJ1+J 

C( 1. JJ)-O.D0 

DO 6 K=1.N 

C( I , JJ)=C( I, JJ)+D( I,K)>^cA(K.JJ) 

RE'rURN 

DO 3 1=1. N 

DO 8 J=1,N 

JJ=NJ1+J 

Hv ! . J; =O.D0 

DO 8 K=1.N 

R( I. J)=R( I.J)+D( I.K)3r:A(K,JJ) 

RETORN 

DO 10 !=1,N 

JJ=NJ1+I 

DO 10 J=l,N 

R( T . J)=O.DO 

DO 10 K=l.N 

R( l.J) = R( I, J)+A(K. JJ)3ecD(K. J) 

RETI.TRN 

END 


F 10 
F 20 

E( IR. IR) . R( IR. IR> F 30 

F 40 
F 60 
F 60 
F 70 
F 80 
F 90 

V 100 
F 1 10 
F 120 
F 130 
F 140 
F 160 
F 160 
F 170 

V 160 
F 190 
F 200 
F 210 
F 220 
F 230 
F 240 
F 250 
F 260 
F 270 
F 280 
F 290 
F 300 
F 310 
F 320 
F 330 
F 340 
F 350 
F 360 
F 370 
F 380 
F 390 
F 400 
F 410 
F 420 
F 430 


SUBROUTINE ZZTT (Z.T3,T4, IR, IC.NB.T5) 

IMPLICIT REAL:«:8(A-H,0-Z) 

DIl'IENSION Z( IR. IG), T3( IR, IR) , T4( IR, IR) , T5( IR. IR) 
COMMON /DIMEN/ N.M 

Mi = NB+N 

no i 1=1, N 

DO I J=l,N 

MM=m+J 

NN=NB+J 

Z( I,NN)=Z( I,MM)-TO( I, J)+T4( I, J) 

T5( I,J)=Z( I,NN) 

RETURN 

END 


G 10 
G 20 
G 30 
G 40 
C 60 
0 60 
G 70 
G 80 
G 90 
G 100 
G 1 10 
C 120 
( 130 

G 14^ 



I 


SUBROUTINE BMULT (C.K.F.NI ,NJ , IR» IC) 

. • COMPUTE PRODUCT OF ’T’ AND PROPER PART OF 
IMPLICIT REALr^t8(A-fl,0-Z) 

DIMENSION X( 1) , F(l), G( IR, IR) 

COMMON /DIMEN/ IMAX, JMAX, MAX 

II=WI/IM\X+1 

JJ=NJ/IM^X+1 

MN=MAX-JJ+1 

NII^CMAX-DsfilMAX-NJ 

NJJ=(PLAX-l):ecIMAX-NI 

CALL AB (C,X,F,NI.NJ, IR, IC) 

IF( n.EO.JJ.AI^D. lI.EO.MM) RETURN 

JF( II.EO. JJ) <X) TO 1 

CALL ABTR ( C,X,F,NJ,NI , IR, lO 

IF( II.EO..MM) RETURN 

CALL AB rG,X,F,NJJ.NII, IR, IC) 

CALL ABTK ( C, X, F, NI I , NJJ, IR, IC) 

RETURN 

END 


*X» TO GET 


H 10 
H 20 
H 30 
H 40 
H 50 
'i 60 
H 70 
ri BO 

H 90 
H 100 
II 110 
H 120 
H 130 
H 140 
H 160 

16«iy 

II 170 
H 180 
H 190 
H 200 


I 

I 


SUBROUTINE AB (T,X,F,NI ,NJ, IR, IC) 
IMPLICIT REAL^8(A-H,0-Z) 
DIMENSION T( IR, IR) , X(l), F( 1) 
COM^ION /DIMEN/ IMAX. JMAX, MAX 


G ... COMPUTES THE PRODUCT OF ’T* AND PROPRE PART OF ’X* AND 
C ADDS TO ’F’ 

DO 1 1=1, IMAX 
II = NI*M 
DO 1 J=l, IMAX 
JJ=NJ+J 

1 F( II)=F( !I)+T( 

RETURN 

END 


10 

20 

30 

40 

50 

60 

70 

80 

90 

100 

110 

120 

130 

140 


C 


1 


SUBROUTINE ABTR (T,X,F,NI ,NJ, IR, IC) 
IMPLICIT REAL:f:8(A-H,0-Z) 

DIMENSION T(IR, IR), X( 1) , FC 1) 
COMMON /DIMEN/ IMAX. JMAX, MAX 


DO 1 1=1. IMAX 

II=NI+I 

DO 1 J=l. IMAX 

JJ=NJ+J 

F( II)=F( II)+T(J, I):#:X(JJ) 

RETURN 

END 


J 10 
J 20 
J 30 
J 40 
J 60 
J 60 
J 70 
J 80 
J 90 
J 100 
J 110 
J 120 


i 

f 

\ 




^IGINAL PArtr 


